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A  multi-block  method  being  explored  is  further  developed  and  investigated 
for  the  simulation  of  Mach  0.8  transonic  turbulent  flow  past  a  wing-fuselage 
configuration  at  -3.0  degree  angle  of  attack.  In  this  method  the  flow  field  of 
interest  is  first  divided  into  six  contiguous  blocks  in  such  a  way  that  each  block  is 
partially  bound  by  a  solid  surface.  Accordingly  the  solid  surface  is  mapped  onto  a 
complete  rectangular  boundary  plane  in  the  computational  domain  for  effective 
calculation  of  algebraic  eddy  viscosity  as  well  as  for  direct  application  of  a  thin- 
layer  Navier-Stokes  code.  Then,  in  the  solution  process,  each  block  is  treated  as 
an  independent  flow  problem,  modeled  by  the  unsteady  Reynolds-averaged  thin- 
layer  Navier-Stokes  equations,  with  the  interface  boundary  conditions  updated  at 


every  time  step.  The  turbulence  closure  model  employed  is  the  Baldwin-Lomax 
eddy  viscosity  model. 

The  multi-block  method  being  investigated  has  a  distinct  advantage  in 
design  application  in  that  a  local  change  to  the  configuration  requires  only  the 
related  block  grid  to  be  regenerated.  However,  the  excessive  distortion  on  block 
domain  transformation,  in  particular  the  wing  blocks,  imposed  by  the  special 
features  of  the  method  makes  it  difficult  to  generate  good  block  grids  for  accurate 
flow  simulation.  Accordingly,  special  measures  and  proper  techniques  for  quaUty 
block  griddings  have  to  be  developed  and  investigated. 

In  this  study  a  one-dimensional  adaptive  grid  generation  technique  based 
on  a  variational  approach  has  been  employed  to  generate  wing  surfaces  adaptive 
to  measured  surface  pressure  gradients.  With  the  solid  surface  and  the 
corresponding  outer  surface  generated,  the  six  block  grids  are  generated  by  an 
algebraic  Two-Boundary  Method  based  on  transfinite  interpolation.  The  blending 
function  used  is  a  Hermite  polynomial  which  can  provide  good  orthogonal  grids 
near  the  surface.  Also,  a  hyperbohc  tangent  clustering  function  has  been  used  to 
provide  effective  grid  point  distribution  and  sufficient  grid  resolution  in  the 
viscous  boundary  layer.  For  the  adaptive  block  grids  generated,  the  Jacobian  of 
transformation  evaluated  by  finite-difference  approximations  is  positive  at  every 
grid  point.  To  gain  further  insight  and  understanding  on  the  multi-block  method 
for  complex  flow  simulation,  five  different  sets  of  six  block  grids  have  been 
generated  and  investigated  for  the  wing-fuselage  flow  simulation.  Numerical 


vi 


results  obtained  show  that  the  multi-block  method  investigated  is  a  very  promising 
approach  which  can  be  further  developed  for  complex  configuration  aerodynamics 
simulation. 
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CHAPTER  I 
INTRODUCTION 


With  the  advent  of  supercomputers  and  the  advancement  in  solution 
algorithms  for  nonlinear  problems,  computational  fluid  dynamics  (CFD)  has 
become  a  powerful  tool  for  analyzing  practical  fluid  flow  problems.  Even  though 
much  progress  has  been  made  in  CFD,  further  advances  are  needed  in  order  to 
effectively  use  it  for  simulating  high  speed  turbulent  flows  past  complex 
aerodynamic  configurations.  A  successful  numerical  simulation  of  flow  past  a 
complex  configuration  requires  not  only  an  accurate  and  robust  flow  solver  but 
also  a  good  grid  system.  In  this  study,  a  computational  method  utilizing  a  multi- 
block  grid  system  is  developed  and  investigated  for  simulating  transonic  turbulent 
flows  past  a  wing-fuselage  configuration. 

We  begin  this  introductory  chapter  with  some  general  background 
information  on  grid  generation.  Afterwards,  a  literature  survey  is  given  for  topics 
relevant  to  the  multi-block  computational  methods  developed  in  this  study.  This 
chapter  concludes  with  a  description  of  the  objectives  of  the  present  study  and  an 
overview  of  this  dissertation. 
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1.1  Background 

Grid  generation  is  concerned  with  the  mapping  of  unequally  distributed 
grid  points  in  the  spatial  domain  onto  a  computational  domain  where  all  of  the 
grid  points  are  uniformly  distributed.  Much  of  the  rapid  progress  in  CFD  can  be 
directly  attributed  to  the  development  of  grid  generation  techniques,  because  they, 
in  conjunction  with  methods  such  as  the  finite-difference  algorithm,  have  enabled 
one  to  study  flow  problems  in  complex  geometries  efficiently  and  accurately  [1-8]. 

Eiseman  and  Erlebacher  [7]  classified  all  possible  grid  systems  that  can  be 
used  by  finite-difference  method  as  follows:  At  the  broadest  level,  a  grid  system 
can  be  classified  as  structured,  unstructured,  or  mixed  depending  upon  how  the 
grid  points  are  connected  to  each  other.  A  structured  grid  system  in  turn  can  be 
classified  as  a  single  grid  or  a  multi-block  grid.  A  single  grid  is  one  that  is  based 
on  a  single  boundary-fitted  coordinate  system  whereas  a  multi-block  grid  is  made 

up  of  two  or  more  single  grids  patched  together  with  each  single  grid  having  a 

different  boundary-fitted  coordinate  system. 

Of  the  grid  systems  mentioned  above,  the  unstructured  and  mixed  grid 

systems  are  more  versatile,  especially  for  compHcated-shaped  spatial  domains; 

but  the  use  of  these  grid  systems  with  a  finite-difference  method  is  still  in  a  state 

of  development  [8,  9].  Presently,  most  finite-difference  methods  use  the 

structured  grid  because  of  its  inherent  compatibility  with  such  methods.  In  this 

research,  the  structured  grid  is  used. 
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Methods  for  generating  grid  systems  can  be  divided  into  two  major  classes- 
differential  equation  methods  and  algebraic  methods.  Differential  equation 
methods  generate  grid  systems  by  solving  a  system  of  partial  differential  equations 
which  describe  how  grid  points  are  to  be  distributed  within  the  spatial  domain. 
Many  of  these  methods  require  a  substantial  amount  of  computational  effort, 
since  the  system  of  partial  differential  equations  that  must  be  solved  is  quasilinear 
and  often  as  complicated  as  the  partial  differential  equations  that  govern  the  fluid 
flow  problem. 

Algebraic  methods  generate  grid  systems  by  interpolating  between 
boundaries  of  the  spatial  domain  and  using  stretching  functions  to  control  the 
distribution  of  grid  points.  Here,  for  the  purpose  of  grid  generation,  interpolation 
means  a  process  of  inferring  grid  point  locations  between  points  where 
information  is  known  exactly.  Since  the  desired  grid  point  locations  are  an 
algebraic  combination  of  known  quantities,  no  partial  differential  equation  needs 
to  be  solved  in  this  grid  generation  process. 

A  "quality"  grid  system,  whether  generated  by  differential  equation  or 
algebraic  methods,  is  usually  one  that  is  reasonably  orthogonal  and  smooth  with 
the  grid  points  concentrated  in  regions  where  they  are  most  needed.  By 
strategically  positioning  the  grid  points,  the  total  number  of  grid  points  needed  to 
obtain  accurate  solutions  can  be  reduced.  This  reduces  computer  memory  and 
CPU  time  requirements. 


In  order  to  cluster  grid  points  where  they  are  most  needed,  it  is  desirable 
to  couple  the  grid  generation  process  with  the  numerical  solution.  Grid 
generation  methods  which  use  knowledge  of  the  solution  to  construct  grid  systems 
are  referred  to  as  solution-adaptive  grid  generation  methods.  In  this  study,  some 
of  the  ideas  utilized  in  generating  solution-adaptive  grids  will  be  employed  in  the 
generation  of  a  quality  multi-block  grid  for  the  wing-fuselage  configuration. 

For  problems  involving  complex-shaped  spatial  domains,  it  is  often  difficult 
to  generate  a  single  grid  system  that  would  be  ideal  for  the  entire  flow  field.  For 
such  problems,  it  is  often  necessary  to  divide  the  spatial  domain  into  blocks  so 
that  an  ideal  grid  system  can  be  effectively  generated  for  each  block.  As 
mentioned  earlier,  a  grid  system  that  is  made  up  of  a  number  of  single  grids 
patched  together  is  known  as  a  multi-block  grid.  In  this  study,  computational 
methods  which  use  multi-block  grids  are  referred  to  as  multi-block  methods. 

1.2  Literature  Survey 

In  keeping  with  the  topics  that  are  relevant  to  the  present  study,  this 
Uterature  survey  is  concerned  with  the  following  three  subject  areas: 

(1)  Algebraic  grid  generation  techniques, 

(2)  Solution-adaptive  grid  generation  techniques,  and 

(3)  Complex  flow  simulations  using  multi-block  grids. 
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1.2.1  Algebraic  Grid  Generation  Techniques 

The  foundation  of  algebraic  grid  generation  can  be  traced  back  to  a  paper 
published  in  1967  by  Coons  [10],  in  which  the  concept  of  blending  functions 
between  curves  and  surfaces  was  introduced.  Gordon  and  Hall  [11]  applied 
Coons'  work  to  numerical  grid  generation,  and  introduced  the  concept  of 
transfinite  interpolation.  Smith  [12]  used  the  findings  of  Coons  as  well  as  Gordon 
and  Hall  to  develop  a  versatile  algebraic  grid  generation  technique  known  as  the 
Two-Boundary  Method.  Eriksson  [13]  and  Vinokur  and  Lombard  [14]  extended 
the  Two-Boundary  Method  so  that  four  boundaries  of  the  spatial  domain  can  be 
mapped  correctly.  The  resulting  method  is  referred  to  as  the  Four-Boundary 
Method.  Yang  and  Shih  [15]  extended  the  Two-Boundary  Method  to  handle 
time-varying  deforming  spatial  domains.  Shih  et  al.  [16]  developed  a  method  for 
generating  parametric  representations  of  boundary  surfaces  needed  in  algebraic 
grid  generation.  They  also  developed  an  efficient  computer  program  to  generate 
single  and  multi-block  grid  systems  for  geometrically  complex  two-  and  three- 
dimensional  spatial  domains  by  using  the  Two-  and  Four-Boundary  Methods, 
1-2.2  Solution- Adaptive  Grid  Generation  Techniques 

Solution-adaptive  grid  generation  techniques  have  been  comprehensively 
surveyed  by  Anderson  [17]  and  Thompson  [2].  Russell  and  Christiansen  [18] 
noted  that  all  solution-adaptive  grid  generation  methods  essentially  attempt  to 
equidistribute  some  weighted  solutions.  Recently,  a  simplified  approach  based  on 
this  equidistribution  concept  is  described  by  Shyy  [19],  in  which  the 


multidimensional  adaptive  procedure  is  replaced  by  a  series  of  one-dimensional 
ones.  This  simplified  approach  reduces  computing  time  considerably. 

The  ideas  of  the  equidistribution  principle  comes  from  a  variation 
formulation  developed  by  Brackbill  and  Saltzman  [20].  They  defined  different 
functionals  to  measure  and  control  different  grid  characteristics  (e.g.,  smoothness, 
orthogonality,  and  clustering)  in  the  spatial  domain.  The  minimization  of  these 
functionals  results  in  a  set  of  partial  differential  equations  which  govern  the  grid 
point  distribution.  These  partial  differential  equations  are  quite  complex,  and 
require  considerable  computational  resources  to  analyze.  Fortunately,  for  many 
flow  problems,  disparate  length  scales  in  the  solution  occur  only  in  one  coordinate 
direction.  For  this  case,  adaptation  is  required  in  only  one  coordinate  direction, 
and  the  partial  differential  equation  reduces  to  an  ordinary  differential  equation. 
The  resulting  ordinary  differential  equation  yields  the  equidistribution  principle. 
Because  of  the  simpUcity  of  the  one-dimensional  variational  formula,  many 
investigations  employ  it  for  higher  dimensions  [6,  19]. 
1.2.3  Complex  Flow  Simulations  Using  Multi-Block  Grids 

Hoist  et  al.  [21]  employed  multi-block  grids  with  the  Euler  equations  and 
thin-layer  Navier-Stokes  equations  to  study  transonic  flows  around  wings.  It  was 
noted  in  their  study  that  proper  interfacing  between  different  blocks  can  be 
difficult  but  is  critical  to  a  successful  simulation.  More  recently,  the  thin-layer 
Navier-Stokes  equations  were  used  vidth  Chimera  overlapped  grids  to  investigate 
subsonic  turbulent  flow  about  the  F-18  fuselage  forebody  and  the  combined  wing- 
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fuselage  [22].  In  their  study,  special  coding  and  interpolation  schemes  were 
required  for  the  overlapped  grid  region,  but  the  computed  results  were  in  good 
agreement  with  flow  visualization  and  surface  pressure  measurements  obtained  in 
flight-test  conditions.  Huband  et  al.  [23]  simulated  the  transonic  turbulent  flow 
field  around  an  F-16A  fighter  configuration  by  solving  the  Reynolds-averaged 
Navier-Stokes  equations  on  a  single  grid.  For  that  single  grid,  orthogonality  of 
grid  lines  at  solid  surfaces  was  not  enforced. 

There  are  a  number  of  multi-block  grid  topologies  proposed  for  the  Euler 
simulation  of  fighter  aircraft  models  and  other  space  vehicle  configurations.  Some 
representative  sample  topologies  can  be  found  in  a  recent  conference  proceeding 
[24].  For  a  transonic  turbulent  flow  simulation,  however,  the  multi-block  grid 
topology  and  the  body-conforming  transformation  for  block  gridding  are  more 
involved  and  require  special  considerations.  This  is  because  boundary  layers  must 
be  properly  discretized  and  a  normal  distance  from  a  grid  point  to  the  solid 
surface  is  required  for  the  calculation  of  eddy  viscosity.  Accordingly,  if  the  flow 
field  is  divided  into  blocks  just  for  ease  in  generating  block  grids  then  special 
measures  will  be  required  for  proper  eddy  viscosity  calculation.  An  improper 
eddy  viscosity  can  be  detrimental  to  accurate  simulation  of  complex  viscous  flow 
phenomena,  even  though  it  may  have  little  effect  on  surface  pressure  distribution. 
Moreover,  a  specific  flow  domain  transformation  is  highly  desirable  for  effective 
use  of  available  grid  points  and  direct  appUcation  of  an  existing  thin-layer  Navier- 
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Stokes  code.  Hence,  a  good  multi-block  grid  topology  for  the  Euler  simulation  in 
general  may  not  be  a  good  one  for  the  turbulent  flow  simulation. 

It  is  clear  that  a  proper  choice  of  grid  topology  and  transformation  is 
critical  to  accurate  and  efficient  simulation  of  turbulent  flow  past  a  complex 
configuration.  Most  recently,  a  novel  multi-block  method  for  complex  turbulent 
flow  simulation  has  been  explored  and  investigated  by  Yang  [25].  In  his  method, 
the  flow  field  was  first  divided  into  a  number  of  contiguous  blocks,  in  such  a  way 
that  each  block  was  partially  bounded  by  a  solid  surface.  Accordingly,  the  soUd 
surface  was  mapped  onto  a  complete  rectangular  boundary  plane  in  the 
computational  space  for  effective  calculation  of  eddy  viscosity  as  well  as  for  direct 
application  of  a  thin-layer  Navier-Stokes  code.  Then,  in  the  solution  process,  each 
block  was  treated  as  an  independent  flow  problem,  modeled  by  the  unsteady  thin- 
layer  Navier-Stokes  equations,  with  the  interface  boundary  conditions  updated  at 
every  time  step.  The  multi-block  method  explored  has  a  distinct  advantage  for 
design  appUcations  in  the  sense  that  a  local  change  to  the  configuration  requires 
only  the  related  block  grid  to  be  regenerated.  This  approach  has  been 
investigated  for  the  simulation  of  Mach  0.8  turbulent  flow  past  a  wing-fuselage 
configuration.  Yang's  study,  however,  was  incomplete  in  that  techniques  were  not 
developed  to  ensure  a  smooth  and  nonsingular  mapping.  It  is  mentioned  in 
passing  that  special  features  of  his  method  has  resulted  in  excessive  distortion  on 
block  domain  transformation  which  makes  it  very  difficult  to  generate  good  multi- 


block  grids  for  flow  simulations.  Despite  these  shortcomings,  Yang's  exploratory 
study  did  demonstrate  the  usefulness  of  the  multi-block  method. 

1.3  Objectives 

As  mentioned  in  the  previous  section,  Yang's  exploratory  study  [25]  showed 
that  the  multi-block  method  has  promise  for  accurate  and  efficient  simulation  of 
high  speed  turbulent  flows  past  complex  configurations.  Therefore,  the  objective 
of  this  study  is  to  further  develop  the  multi-block  method  for  complex  ! 
configuration  flow  simulations.  More  specifically,  the  objectives  are  as  follows. 
First,  develop  an  efficient  multi-block  computational  method  for  the  simulation  of 
transonic  turbulent  flow  past  a  wing-fuselage  configuration  that  ensures  a  smooth 
and  nonsingular  mapping.  Second,  assess  the  importance  of  grid  resolution  on 
solution  accuracy.  And  finally  third,  study  the  physics  of  simulated  flow  problem. 

1.4  Overview 

In  the  next  chapter,  the  flow  problem  chosen  for  study  and  the  governing 
equations  suitable  for  solving  the  problem  are  described.  In  Chapter  III,  the  grid 
topology  and  flow  solver  are  described.  Chapter  IV  presents  the  multi-block  grid 
generation  method.  Solutions  obtained  by  using  the  method  developed  in  this 
study  are  presented  in  Chapter  V.  Finally,  Chapter  VI  provides  a  brief  summary 
of  the  major  points  of  this  dissertation. 


CHAPTER  n 
FORMULATION  OF  PROBLEM 


This  chapter  describes  the  flow  problem  to  be  studied  and  its  governing 
equations.  The  flow  problem  is  described  in  Section  2.1,  and  the  equations  which 
govern  the  problem  are  described  in  the  following  order.  In  Section  2.2,  the 
nondimensionalized  Navier-Stokes  equations  are  given.  In  Section  2.3,  these 
equations  are  transformed  to  generalized  coordinates  to  facilitate  the  numerical 
method  of  solution.  In  Section  2.4,  the  governing  equations  in  generalized 
coordinates  are  simplified  to  a  thin-layer  form  for  reducing  computational  cost. 
In  Section  2.5,  the  turbulence  model  used  is  described  in  detail.  In  Section  2.6, 
boundary  conditions  are  specified.  Finally,  initial  conditions  are  described  in 
Section  2.7. 

2.1  Description  of  Problem 

The  wing-fuselage  configuration  shown  in  Figure  2.1  is  selected  to  further 
develop  the  multi-block  computational  method.  The  flow  problem  considered  is  a 
transonic  turbulent  flow  of  Mach  number  (M  J  0.8  past  the  configuration  at  an 
angle  of  attack  (a)  of  -3°.  According  to  the  wind-tunnel  test  conditions,  the  free- 
stream  stagnation  temperature  is  609°R,  and  the  Reynolds  number  (Re)  is 
2.2x10^. 
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By  taking  the  diameter  of  the  largest  cross-section  of  the  fuselage  (D)  as  a 
reference  length  scale,  the  length  of  fuselage  is  about  7D;  the  wing  span  is  about 
3. 3D;  the  chord  length  at  the  wing  root  is  about  1.4D.  Additional  details  of  the 
wing-fuselage  geometry  can  be  found  in  the  work  of  Yang  [25]  and  will  not  be 
repeated  here. 

Here,  it  is  assumed  that  the  flow  field  is  symmetric  with  respect  to  the 
plane  located  at  z  =  0.  With  this  assumption,  the  spatial  domain  of  interest  is  the 
region  shown  in  Figure  2.2.  The  distance  between  the  inflow  and  outflow 
boundaries  is  37D,  and  the  distance  between  the  plane  of  symmetry  at  the 
fuselage  and  the  outer  boundary  is  about  lOD. 

The  wing-fuselage  problem  just  described  is  selected  for  the  present  study 
primarily  because  of  the  following  reasons.  First,  experimental  data  of  the 
pressure  coefficient  on  the  wing  are  available  for  evaluating  the  computational 
method  developed  and  investigated  in  this  study.  Second,  this  is  a  stringent  test 
problem  since  the  turbulent  flow  field  is  complex,  involving  shock  wave-boundary 
layer  interaction  and  possible  flow  separation.  Third,  this  problem  is  of  strong 
current  interest  since  the  flow  configuration  models  the  next  generation  transonic 
propfan  transport  aircraft.  Finally,  the  geometry  of  the  wing-fuselage 
configuration  is  complex  enough  to  test  the  grid  generation  techniques  of  the 
computational  method  being  developed. 
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Figure  2.2  Spatial  domain 
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2.2  Nondimensionalized  Navier-Stokes  Equations 


The  equations  which  govern  the  flow  problem  described  in  the  previous 
section  are  the  Navier-Stokes  equations  for  ideal  gases  without  either  body  forces 
or  external  heat  addition.  With  a  superscript  asterisk  denoting  a  dimensional 
variable,  these  equations  can  be  nondimensionalized  by  dividing  all  lengths  by  a 
characteristic  length  scale  (L*),  all  velocities  by  the  free-stream  speed  of  sound 
(a„*),  density  by  the  free-stream  density  (Po.*),  time  by  L'/slJ  ,  internal  energy 
per  unit  mass  by  a„'^  pressure  by  p„*a„*^,  and  dynamic  viscosity  by  The 
nondimensionaUzed  form  of  these  equations  can  be  written  in  the  following  strong 
conservation-law  form  [26] 

— -±  + — + — + —  =  Re  (  -+  -+  ! 


Re-\. 

dt     dx    dy    dz  dx.     dy  dz 


) 


where  the  Reynolds  number  is  defined  as  Re  =  p.'a„'L*/ki„*  ,  and 


Q  = 


p 

pu 
pv 
pw 
e 


(2.1) 


E  = 


pu 

pv 

pw 

pu^+p 

pvu 

pwu 

puv 

F  = 

pv^+p 

G  = 

pwv 

puw 

pvw 

pw^+p 

(e+p)u. 

.(e  +p)v 

(e+p)w 
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xz 

K 

K 
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The  elements  of  vector  Q  are  the  nondimensionalized  dependent  variables.  In  the 
above  equations,  p  is  nondimensionalized  density;  u,  v,  and  w  are,  respectively, 
the  X-,  y-,  and  z-  components  of  nondimensionalized  velocity;  and  e  is  the 
nondimensionalized  total  energy  per  unit  volume  defined  by 

e  =  p[ei  +  ^(u^+v^+w^)] 

The  internal  energy    appearing  in  the  above  equation  is  given  by 

=  T 
~  Y(Y-l) 

where  the  ratio  of  specific  heats  y  =  C^'/Cj  =  1.4  for  air,  with  Cp*  and  Cj  are 
the  specific  heats  at  constant  pressure  and  volume,  respectively.  The 
nondimensionalized  pressure  p  and  temperature  T  which  appear  in  the  governing 
equations  can  be  expressed  in  terms  of  the  dependent  variables  by 

p  =  (Y-l)[e--p(u2+v2+w2)] 
T  =  y(y-1)[^-V+v2-.w2)] 
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The  components  of  viscous  stress  tensor  are  given  by 

^xx  =  ^(Ux+Vy+w^)+2jiu^  =  |ii(2u,-Vy-wJ 

^zz  =  ^(Ux-^Vy+wJ+2iiw^  =  |p(2w^-u^-vp 


yx 


T      =  11(11  +V  )  =  T 

T      =  u(v  +W  )  =  T 
yz       f^V'z       y/  zy 

T      =  |J.(w  +U  )  =  T 

zx       f^v    X      z/  xz 

where  \l  is  the  nondimensionahzed  dynamic  viscosity,  and  from  Stokes'  hypothesis 
the  second  coefficient  of  viscosity  A  is  set  equal  to  -2/3  \i. 

The  three  components  of  p  which  represent  the  diffusion  of  mechanical 
and  thermal  energy  are  given  by 

B    =  UT    +VT    +WT   I 

=  UT    +VT    +WT  *   

y  >  Pray 

where  the  Prandtl  number  Pr  =  CpV/k'  =  0.72  for  air  at  standard  conditions, 
and  k'  is  the  thermal  conductivity. 
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2.3  Transformed  Navier-Stokes  Equations 


To  enhance  the  accuracy  and  efficiency  of  a  numerical  scheme  and  to 

simphfy  implementation  of  boundary  conditions  for  complex  geometry,  Eq.  (2.1)  is 

transformed  from  the  x-y-z-t  Cartesian  coordinates  to  general  curvilinear 

coordinates  ^-ri-C-f  by  the  following  relations 

$  =  $(x,y,z,t) 

n  =  n(x,y,z,t) 
C  =  C(x,y,z,t) 
T  =  t 


The  resulting  transformed  nondimensionalized  equations  which  maintain  the 
strong  conservation  law  form  are  given  by  Anderson  et  al.  [26] 


where 


dx   di  dx]  dc         ^di    dx]  dz 


Q  = 


p 

pU 

pu 

puU  +  ^jp 

pv 

t  =  J-i 

pvU  +  $yp 

pw 

pwU  +  $  J) 

e 

(e+p)U-$j) 

(2.2) 


F  =  J-i 


pV 

pW 

puV  +  Tl^p 

puW  +  f^p 

pVV+llyP 

G  =  J-i 

PVW   +  CyP 

PWV  +  TIJ) 

pwW  +  C  J) 

(e+p)V-tij) 

(e+p)W-C,p 
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A  1 


0 

£  T     +£  T     +£  T 
*x  XX     *y  xy     *z  xz 

^x  yx    *y  yy  yz 

£  T    +£  T    +£  T 

'x  zx       y  zy     'z  zz 

^x^x    ^y^y  '''  ^z^z 


0 


nx'^xx-^nyi:xy-^^z'>^xz 

nxV^^y^yy^'^zV 

nx'^zx-^Vzy'-^z'^zz 

^xPx-^nyPy+nA. 


0 

'X  XX       y  xy  ^z  xz 

„      J     ^x  yx    *y  yy  ^z  yz 

*x  zx     ^y  zy  ^z  zz 

^x^x  *  ^yPy  ^z^z 


In  the  above  equations,  U,  V,  and  W  are  components  of  the  contravariant  velocity 
normal  to  constant  £,  r\,  and  C  planes,  respectively,  and  are  given  by 

U  =  $t*^xU  +  5yV  +  $,w 

V  =  tit+n^u+tiyV+n^w 

W  =  Ct  +  ^xU  +  CyV  +  C^w 
The  inverse  of  the  Jacobian  is  given  by 

j-i  ^  d(x,y,z) 
a(£,n,C) 

=  ^«(y.zc-yc^)  -^(y^Zc-y(Z^)  ^x^Cy^z^-y^z^) 
The  metric  coefficients  which  appear  in  the  above  equations  are  given  as 

=  J(y.zc-z^y^)  =  J(y^Zj-Zjyj) 

£y  =  J(z^x^-x^Zj)  tiy  =  J(z^x,-x^Zj) 

5z  =  J(V(-yr,^c)        %  =  J(xcy£-y(Xj) 
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=  HhYr^-yi^)         =  -x.C,-y,Cy-z,C, 
It  is  noted  that  the  Cartesian  derivatives  for  the  components  of  the 
nondimensionalized  stress  tensor  and  heat  flux  vector  have  to  be  expanded  in 
terms  of  i,  x],  and  C  via  chain-rule  relations  such  as 

y     ^y  t    'y  1   ^y  C 

For  instance, 

=  Ji(Uy+v,) 

=  ji(uj$y+u^ny+UjCy+Vj$^+v^nx^VjO 

2.4  Thin-Layer  Navier-Stokes  Equations 

For  the  high  Reynolds  number  transonic  flow  considered  in  this  study,  the 
viscous  effects  are  mostly  confined  to  a  thin  layer  near  solid  surfaces  where  the 
dependent  variables  undergo  rapid  gradient  changes.  Accordingly,  only  viscous 
terms  in  the  direction  normal  to  the  surface  (the  C  direction  in  the  present 
applications)  need  to  be  retained.  This  leads  to  the  thin-layer  approximation 
which  neglects  all  viscous  derivatives  in  the  directions  parallel  to  the  solid  surface 
(the  $  and  n  directions)  and  maps  the  soUd  surface  onto  a  constant  C  plane.  The 
thin-layer  approximation  retains  all  three  of  the  momentum  equations  and  makes 
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no  assumption  about  the  pressure  [26].  Thus,  one  principal  advantage  of  the  thin- 
layer  approximation  is  the  retaining  of  the  normal  momentum  equation  so  that 
pressure  can  vary  through  the  boundary  layer. 

After  removing  the  notation    in  Eq.  (2.2),  the  thin-layer  Navier-Stokes 
equations  in  general  curvilinear  coordinates  S-n-C-t  can  be  written  as  [26] 


where  viscous  derivatives  in  the  direction  normal  to  the  solid  surface  (the  C 
direction)  are  retained  and  collected  into  the  vector  S  given  by 

0 


ar  a$  an  ac  ac 


(2.3) 


s  =  j-^ 


m^Uj+m^C, 
m^Vj+nijCy 
m,w^+m2C, 
m^nig  +  m2(C^u  +  CyV  +  C^w) 


where 


^2  =  |(fxU(  +  CyV^  +  C,Wj) 

m3  =  ^(u2+v2+w2)^  +  Pr-\y-l)-V)c 


2.5  Turbulence  Model 


Since  present  computing  capabilities  prohibit  direct  or  large-eddy 
simulations  of  turbulent  flows  over  a  complex  configuration,  the  thin-layer  Navier- 
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Stokes  equations  are  Reynolds-averaged  to  reduce  computational  cost.  This 
involves  decomposing  the  dependent  variables  in  Eq.  (2.3)  into  mean  and 
fluctuating  components,  and  then  ensemble  averaging  the  entire  equation.  In  the 
averaging  process,  new  unknowns  are  introduced  and  must  be  closed  by  a 
turbulence  model.  The  two-layer  algebraic  eddy  viscosity  model  of  Baldwin  and 
Lomax  [27]  is  considered  in  this  study. 

By  utilizing  the  approach  just  described  for  modelling  turbulence  and 
neglecting  density  fluctuation  effects,  the  Reynolds-averaged  thin-layer  Navier- 
Stokes  equations,  cast  in  strong  conservation  law  form  with  generalized 
coordinates,  have  the  same  form  as  Eq.  (2.3).  However,  all  dependent  variables 
should  now  be  interpreted  as  Reynolds-averaged  mean  quantities.  Also,  the 
coefficient  of  viscosity  \i  is  replaced  with  n  +  Hx  ,  and  the  thermal  conductivity  k  is 
replaced  with  k+k^  ,  where      and  k^  are  the  turbulent  eddy  viscosity  and 
turbulent  thermal  conductivity,  respectively.  The  turbulent  thermal  conductivity 
can  be  expressed  in  terms  of  the  eddy  viscosity  using  the  turbulent  Prandtl 
number  Pr^. ,  which  is  defined  as  Ftj  =  ji^'Cp'/kT*  and  is  taken  as  0.9  here. 
Below,  the  method  used  to  model  jij  is  described.  In  this  description  and  for  the 
rest  of  this  section,  the  superscript  *  which  is  used  to  designate  the  dimensional 
quantities  will  be  dropped. 

As  proposed  by  Baldwin  and  Lomax,  the  eddy  viscosity  coefficient  is 
modelled  in  a  two-layer  fashion  as  shown  below 

(^''T)inner  '  ycrossovcr 
(^^T)outer  '  y ycrossovcr 
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where  y  is  the  normal  distance  from  the  wall,  and  ycrossover  is  the  smallest  value  of 

y  at  which  (llT)in„er  =  (Ht)  outer  • 

The  eddy  viscosity  coefficient  in  the  inner  region  is  based  on  the  Prandtl- 
van  Driest  formulation  [28]  and  is  expressed  as 

(Minner  =  PI^M 

where  |w|  is  the  magnitude  of  vorticity  given  by 

|w|  =       -  v,)^  +  (v,  -  Wy)2  +  (w^  -uj2 

The  parameter  ?  is  the  mixing  length  corrected  with  the  van  Driest  damping 
factor  [29],  that  is 

f  =  ry[l-exp(^)] 
A* 

In  the  above  equation,  the  von  Karman  constant  k"^  =  0.4;  the  damping  constant 
A*  =  26;  and  y*  is  law-of-the-wall  coordinate  defined  as 

in  which  ,  ,  and  ji^  are  the  local  shear  stress,  density,  and  laminar  viscosity 
evaluated  at  the  wall,  respectively. 

The  eddy  viscosity  coefficient  in  the  outer  region  follows  the  Clauser 
formulation  [30]  as  shown  below 

(t^Auter  =  ^^CcppF^^Fj^g0(y) 
where  the  Clauser  constant  K  =  0.0168,  and  Ccp  =  1.6.  The  parameter  is 
found  via 
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^WAKE  ^  ymax^max  '      QvKymax^DIF/Fmax  } 

In  the  above  equation,        =  0.25,  and  the  quantities       and  are 
determined  from  the  function 

F(y)  =y|w|[l-exp(^)] 

A* 

where  F^^^  is  the  maximum  value  of  F(y)  and  y^^^  is  the  value  of  y  at  which  F^^^ 
occurs.  The  exponential  van  Driest  damping  term  exp(-y^/A*)  in  the  above 
equation  is  set  equal  to  zero  in  wakes.  The  difference  between  maximum  and 
minimum  total  velocity  within  the  profile  at  a  fixed  $  and  r|  location  is 

Ui^ip  =  (\/u^.v^.w2  L-  (/u^.v^.w^  L 

The  second  term  in        is  taken  to  be  zero  except  in  wakes.  The  function 
FklebCy)  is  the  Klebanoff  intermittency  factor  given  by 

FKLEB(y)  =  [l-5.5(52^^)'*]-i       and       C^eb  =  0.3 

ymax 

The  Baldwin  and  Lomax  model  has  the  advantage  of  using  the  distribution 
of  vorticity  to  determine  the  length  scales  for  both  the  inner  and  outer  layers, 
thereby  avoiding  the  necessity  of  finding  the  outer  edge  of  the  boundary  layer. 
This  model  has  been  used  extensively  to  simulate  a  variety  of  transonic  turbulent 
flow  problems,  including  the  satisfactory  results  for  projectile  flow  problems 
[31-33].  This  is  the  prime  reason  why  this  model  is  employed  in  the  present  study. 
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2.6  Boundary  Conditions 

In  order  to  obtain  solutions  to  the  Reynolds-averaged  thin-layer  Navier- 
Stokes  equations,  the  values  of  certain  dependent  variables  p,  u,  v,  w,  and  e  must 
be  specified  on  some  portion  of  the  boundaries.  The  boundary  conditions  (BCs) 
imposed  are  as  follows.  At  all  solid  walls,  the  no-slip  condition  is  applied,  and  the 
first  derivative  of  density  normal  to  the  wall  is  set  to  zero.  At  the  inflow,  outflow, 
and  far-field  boundaries,  characteristics  entering  the  domain  determine  the 
number  of  permissible  BCs  and  appropriate  Riemann  invariants  should  be 
imposed  as  BCs.  For  the  present  problem,  the  following  simpUfied  BCs  are  used 
and  found  to  be  adequate.  At  the  inflow  and  far-field  boundaries,  all  flow  >, 
variables  are  kept  constant  at  the  free-stream  conditions.  Since  the  flow  is 
subsonic  at  these  boundaries,  this  set  of  BCs  overly  constrained  the  problem.  This 
overspecification  is  compensated  by  extrapolating  every  flow  variable  at  the 
outflow  boundary.  Finally,  the  boundary  conditions  at  the  symmetry  plane  are 
imposed  by  using  reflected  image  grid  points  beyond  the  symmetry  boundary. 

The  boundary  conditions  described  above  are  those  permitted  by  the 
governing  equations.  In  practice,  additional  boundary  conditions  (known  as 
numerical  boundary  conditions)  are  needed  by  the  numerical  method  of  solution. 
The  required  numerical  boundary  conditions  along  with  details  of  the  boundary 
conditions  outhned  above  are  described  in  Section  3.2.2. 
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2.7  Initial  Conditions 


The  initial  conditions  employed  in  this  study  are  as  follows.  All  dependent 
variables  are  set  equal  to  their  corresponding  free-stream  values  at  every  point 
within  the  spatial  domain.  In  order  to  ease  the  discontinuity  of  the  dependent 
variables  on  the  solid  surface  due  to  the  no-slip  boundary  condition  imposed 
there,  the  boimdary  condition  is  turned  on  slowly  over  a  period  of  time.  Thus,  in 
the  thin-layer  Navier-Stokes  flow  solver,  the  boundary  conditions  for  every 
dependent  variable  at  the  solid  surface  is  scaled  by  slowly  increasing  a  polynomial 
function  SL  such  that 

C       q  =  SLxq,,-H(l-SL)q.  ( 

in  which 

SL  =  (10-15x6  +  6x62)6^ 

0  _  (NC-1)  -  '  - 

30 

In  the  above  equations,  q^,.  is  the  boundary  conditions  with  the  no-slip  condition 
imposed,  and  q.  is  the  boundary  conditions  based  on  uniform  free-stream 
conditions.  The  parameter  NC  is  the  time  step,  and  NC  =  1,  2,  When 
NC^30,q  =  q,,. 


CHAPTER  III 
NUMERICAL  METHOD  OF  SOLUTION 


In  this  study,  a  multi-block  computational  method  is  developed  and 
investigated  for  studying  the  transonic  turbulent  flow  problem  formulated  in 
Chapter  II.  The  multi-block  computational  method  being  developed  involves 
three  components.  These  components  are  grid  topology,  grid  generation,  and  flow 
solver.  The  first  component,  grid  topology,  is  concerned  with  partitioning  of  the 
spatial  domain  into  contiguous  blocks.  The  goal  here  is  to  create  a  proper  grid 
topology  suitable  for  the  wing-fuselage  configuration  so  that  quality  grid  systems 
can  be  generated  for  turbulent  flow  simulations.  The  second  component,  grid 
generation,  is  concerned  with  the  generation  of  grid  systems  for  each  block  of  the 
partitioned  domain.  The  key  feature  here  is  that  the  grid  generation  process 
takes  into  account  the  nature  of  the  flow  field.  The  third  and  final  component, 
flow  solver,  is  concerned  with  obtaining  numerical  solutions  to  the  governing 
equations  formulated  in  Chapter  II. 

In  this  chapter,  we  first  present  an  overview  of  the  method  of  solution. 
Afterwards,  a  detailed  description  is  given  of  the  grid  topology  and  flow  solver. 
The  details  of  the  grid  generation  process  are  presented  in  the  next  chapter. 
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3.1  Grid  Topology 


In  numerical  simulations,  the  use  of  a  single-block  grid  for  simple 
geometries  can  lead  to  efficient  solution  procedures.  However,  for  flow  past  a 
complex  three-dimensional  configuration,  such  as  a  wing-fuselage  configuration 
with  sharp  comers  and  deep  cavities,  the  generation  of  an  acceptable  single-block 
grid  for  effective  turbulent  flow  simulation  is  generally  not  feasible.  One 
alternative  is  to  construct  several  single-grids  for  different  parts  of  the  spatial 
domain  and  then  patch  them  together.  Grids  thus  constructed  are  known  as 
multi-block  grids. 

Here  it  is  noted  that  the  topology  of  a  multi-block  grid  for  turbulent  flows 
past  complex  configurations  is  more  involved  than  those  for  inviscid  or  laminar 
flows.  This  is  because  turbulence  models  require  information  regarding  the 
distance  between  points  in  the  interior  of  the  domain  and  the  solid  surface.  Thus, 
solid  surfaces  of  each  block  must  be  mapped  onto  a  boundary  plane  in  the 
computational  space  for  ease  in  calculating  the  distance  as  well  as  for  direct 
application  of  the  thin-layer  Navier-Stokes  equations. 

A  six-block  grid  topology  which  decomposed  a  single  block  H-grid  provided 
by  the  NASA  Ames  Research  Center  for  the  wing-fuselage  configuration,  was 
developed  and  described  in  a  previous  study  [25].  That  grid  topology  is  adopted 
here.  As  shown  in  Figure  3.1,  the  wing  is  divided  into  2  blocks,  and  the  rest  of 
the  flow  field  is  divided  into  4  additional  blocks. 
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(b)  Cross-sectional  view 
Figure  3.1  A  six-block  composite  grid  topology 
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By  doing  so,  the  flow  field  is  divided  into  six  contiguous  blocks  with  each  block 
partially  bounded  by  a  solid  surface.  The  transformation  between  physical  space 
and  computational  space  of  irmer  and  outer  surfaces  for  Block  1  to  6  is  shown  in 
Figure  3.2.  The  solid  surface  is  mapped  completely  onto  a  rectangular  boundary 
plane  C  =  0  in  the  computational  space. 

The  selected  grid  topology  allows  for  easy  application  of  the  thin-layer 
Navier-Stokes  equations  as  well  as  for  effective  calculation  of  the  algebraic  eddy 
viscosity.  However,  this  grid  topology  does  not  allow  easy  generation  of  quality 
grids  for  each  block.  For  instance,  for  the  nose  block.  Block  1,  shown  in  Figure 
3.2(a),  the  boundary  plane  AEJF  in  the  computational  domain  degenerates  into  a 
single  line  in  physical  space  and  care  must  be  taken  when  dealing  with  the 
boundary  conditions,  as  will  be  detailed  later.  Also,  for  the  wing  blocks.  Block  3 
and  Block  4,  shown  in  Figure  3.2(c,d),  four  of  the  six  bounding  surfaces  in  the 
computational  domain,  CLQR,  LNTQ,  RSTQ,  and  RSHC,  form  a  nearly  flat 
surface  in  the  physical  space.  Therefore,  special  measures  and  care  in  grid 
generation  techniques  are  required  to  overcome  the  difficulty  resulting  from  the 
drastic  change  in  the  domain  transformation. 

Another  difficulty  which  must  be  dealt  with  in  the  generation  of  multi- 
block  grids  is  the  geometric  complexity  of  the  physical  boundary  surfaces.  On  a 
multi-block  grid,  the  overall  domain  is  subdivided  and  each  individual  grid  must 
be  generated  subject  to  specified  boundary  constraints.  As  shown  in  Figure  3.3, 
the  partitioning  results  in  common  boundaries  that  link  the  information  between 


(b)  Block  2 

Figure  3.2  Block  transformation  between  physical  space  and  computational  space 
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(d)  Block  4 
Figure  3.2  (continued) 
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BLOCK  6 


(f)  Block  6 
Figure  3.2  (continued) 


(b)  Outer  surfaces 
Figure  3.3  The  six-block  surface  interfaces 


neighboring  blocks.  The  interface  boundary  information  must  be  supplied 
without  degrading  numerical  accuracy  and  convergence. 

In  this  study,  a  method  is  developed  to  overcome  the  above  described 
shortcomings  associated  with  the  block  topology  selected.  The  details  of  this 
method  are  presented  in  Chapter  IV. 

Before  leaving  this  section,  we  note  two  things.  First,  in  the  solution 
process,  each  block  is  solved  one  at  a  time,  one  after  the  other.  Each  block  is 
treated  as  an  independent  flow  problem  and  governed  by  the  same  thin-layer 
Navier-Stokes  equations  with  the  interface  boundary  conditions  updated  by  a 
distance  weighted  interpolation  at  the  end  of  each  time  step.  Details  of  such  an 
interpolation  procedure  will  also  be  given  later.  Second,  as  will  be  described  in 
Chapter  IV,  the  computational  domain  in  the  ^-ti-C-t  coordinate  system  is 
replaced  by  a  system  of  equally  spaced  grid  points  and  the  indices  J,  K,  and  L 
denote  the  location  of  each  grid  point.  The  J,  K,  and  L  indices  correspond  to 
x],  and  C  directions,  respectively.  Here,  $  is  along  the  streamwise  direction,  x]  is 
along  cross-section  direction,  and  C  is  normal  to  solid  surface. 

3.2  Flow  Solver 

As  noted  earlier,  the  flow  solver  is  concerned  with  obtaining  numerical 
solutions  to  the  governing  equations  formulated  in  Chapter  11.  In  this  section,  the 
algorithm  employed  is  described  in  detail  along  with  the  required  numerical 
boundary  conditions. 
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3.2.1  TVD  Algorithm  f  ■ 

Solutions  to  the  governing  equations  formulated  in  Chapter  II  can  only  be 
obtained  by  numerical  methods.  Successful  numerical  simulations  of 
aerodynamics  problems  depend  strongly  upon  an  accurate  and  robust  algorithm 
and  a  high  quality  grid  system. 

A  versatile  computer  program  known  as  ARC3D  was  developed  by  NASA 
Ames  Research  Center  to  analyze  steady  or  unsteady,  inviscid  or  viscous  three- 
dimensional  high  speed  compressible  flows  past  arbitrarily-shaped  geometric 
objects  [34].  An  original  version  of  ARC3D  approximates  the  thin-layer  Navier- 
Stokes  equations  for  ideal  gas  in  strong  conservation-law  form  by  using  an  ADI- 
like  algorithm  developed  by  Beam  and  Warming  [35].  In  the  Beam  and  Warming 
scheme,  central  differencing  formulas  are  used  for  both  convection  and  diffusion 
terms.   As  a  result,  second-order  implicit  and  fourth-order  explicit  artificial 
dissipations  terms  are  needed  to  maintain  numerical  stabiUty. 

In  a  previous  research  effort  to  develop  solution-adaptive  computational 
methods  for  transonic  turbulent  flow  past  a  secant-ogive-cylinder-boattail 
projectile  [31,  32],  an  axisymmetric  version  of  ARC3D  was  used  as  the  base  code 
for  numerical  simulation.  In  that  effort,  it  was  found  that  for  complex  flows  with 
shock  wave-boundary  layer  interactions  and  large  separated  regions,  the  Beam 
and  Warming  scheme  is  highly  sensitive  to  the  grid  system  and  has  convergence   '  - 
problems.  In  order  to  overcome  these  difficulties,  the  solution  algorithm  of  the 
axisymmetric  Navier-Stokes  code  was  modified  to  an  implicit  Total  Variation 
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Diminishing  (TVD)  scheme  developed  by  Yee  [36,  37].  The  TVD  algorithm 
employed  is  very  similar  to  the  Beam  and  Warming  scheme  except  that  more 
sophisticated  dissipation  terms  are  used.  As  reported  in  [32],  the  axisymmetric 
TVD  algorithm  is  indeed  more  robust  than  the  Beam  and  Warming  scheme. 
Hence,  an  ARC3D  code  had  been  modified  to  a  TVD  algorithm  and  numerical 
experiments  conducted  showed  that  the  TVD  scheme  in  general  gives  more 
accurate  solutions  for  transonic  turbulent  flows  past  a  projectile  model  at  10° 
angle  of  attack  [33].  Accordingly,  the  modified  ARC3D  code  with  the  TVD 
algorithm  is  employed,  and  the  details  of  this  algorithm  are  described  below  in 
the  following  steps: 

(1)  Temporal  differencing  and  linearization; 

(2)  Approximate  factorization; 

(3)  Spatial  differencing  and  artificial  dissipation; 

(4)  Operator  splitting  and  solution  algorithm. 
Temporal  differencing  and  linearization 

In  the  TVD  thin-layer  Navier-Stokes  code,  the  time  derivative  of  Eq.  (2.3) 
is  approximated  by  a  first-order  accurate  Euler  implicit  formula 

AQ°  =  Q""^-Q°  ^  ■ 

=  At(^)     +0[(At)2]  ^, 

Here,  Q"  =  Q(nAx),  and  Ax  is  the  time  step  size.  In  the  above  equation,  Q"  is 
assumed  to  be  known  and  Q""'Ms  the  unknown  sought.  By  substituting  Eq.  (2.3) 
into  Eq.  (3.1),  and  using  delta  formulation,  we  obtain 
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AQ-  =  -AT[|^.|?.^-Re-H|)]-U0[(Ax)2]  (3.2) 
d$    dx]    dC  dC 

In  the  above  equation,  the  flux  vectors  E,  F  and  G  at  time  level  n  + 1  are 
nonlinear  functions  of  the  unknown  vector  Q""^^  This  nonlinear  nature  can  be 
removed  by  using  Taylor  series  expansion  about  Q",  as  shown  below 

E"*^  =  E"+A"AQ"+0[(At)2] 

=  P  +  B"AQ"+0[(At)2]  (3.3) 

G""^  =  G"  +  C"AQ"+0[(At)2] 

where  A,  B,  C  are  the  inviscid  flux  Jacobian  matrices  defined  by 

A  =  ^,       B  =  ^,       C  =  ^ 
dQ  dQ  dQ 

The  viscous  flux  vector  S  at  time  level  n+1  which  involves  gradients  of  the 

dependent  variable  and  the  Jacobian,  is  linearized  as  suggested  by  Steger  [38]  with 

the  assumption  that  viscosity  coefficient  \l  and  thermal  conductivity  k  are  locally 

independent  of  the  unknown  vector  Q"+\  so  that  elements  of  S  are  expressed  as 

s.  =  j-X^ 

where  is  independent  of  Q"*^  and  pj  is  a  function  of  Q'''^^  These  elements  are 
then  linearized  by  Taylor  series  expansion  as 

Sr  =  S.%J-«r|(j:^"Aq^).0[(Ax)^] 

The  resulting  expression  for  S"^^  can  be  written  in  vector  form  as 

gn.i  ^  s"+J-^M"JAQ"  +  0[(At)2]  (3.4) 


37 

where  M  is  the  viscous  flux  Jacobian  matrix  with  elements  given  by 

After  substituting  Eqs.  (3.3)  and  (3.4)  into  Eq.  (3,2),  we  obtain 
{1  + At[^+^ +^ -Re-^(^llM^)]}"AQ» 

di  dt)  dc  dc 

=  -Ax[f  .^.f -Re-^(|)r.O[(Ax)^] 
d$    dT|    dC  dC 


(3.5) 


where  I  is  a  five  by  five  identity  matrix. 
Approximate  factorization 

Solving  Eq.  (3.5)  directly  for  the  unknown  AQ"  will  involve  the  inversion  of 
a  very  large  system  of  coupled  equations  which  is  rather  tedious  and  expensive. 
Beam  and  Warming  [35]  employed  an  approximate  factorization  procedure  to 
replace  the  original  system  by  a  number  of  much  smaller  systems  of  equations. 
By  following  Beam  and  Warming  [35],  Eq.  (3.5)  can  be  factored  as 
[I  +  A T§]"[I  +  ^x^r^  At[^  - Re-i(^J''MJ 

(3  6) 

The  error  introduced  by  this  factorization  is  second-order  in  time. 
Spatial  differencing  and  artificial  dissipations 

The  spatial  derivatives  appearing  in  Eq.  (3.6)  have  to  be  replaced  by  finite- 
difference  formulas.  All  diffusion  terms  are  approximated  by  second-order 
accurate  central  differencing  formulas.  All  convection  terms  are  approximated  by 
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using  a  TVD  type  differencing.  The  resulting  finite-difference  equations  have  the 
following  form; 

[I  +  At6jA  +  DjJ"[I  +  At6  B+D.  ]"{!  + AT[6,C-Re-^fi Jr^MJ)] +D„}"AQ" 

(3.7) 

=  -AT[6jE  +  6^F  +  6jG-Re-^6^S]"  +  Dg 
where  6  denotes  second-order  central  differencing.  The  implicit  operators  and 
explicit  dissipation  term  appear  in  the  above  equation  are  given  by 

2'^  2 

2  2 

Die  =  -0-5At(Q,,,,,i-Q,^,.i)  (3.8) 
De  =  -0.5At(T^^i$^^i-T^_iO^_,.T^^i$^^, 

2        2  2        2  2  2 


where 


"j.i,K,L  =  (diag[-maxT(Aj')]}^  ,A^^i 

2  2  2 

"j,K.i,L  =  {diag[-maxT(X™)]}^  lA^^i  (3.9) 

2  2  2 

"j.K,L.2  =  {diag[-maxT(A.J)]}^  ,A^^i 

2  2  2 

Here,  diag(z)  denotes  a  diagonal  matrix  with  diagonal  elements  z,  and  i|r(z)  is  an 
entropy  correction  function  given  by 


|z|  ,|z|^e=0.25 

T(z)=1,  2    2,  (3-^^) 
^^-pj.  ,|z|<€=0.25 


39 

and  Xj"*,  A^™,  Xj"  are,  respectively,  the  m"'  eigenvalue  of  the  Jacobian  matrices 

A,  B  and  C  of  Eq.  (3.3).  In  Eq.  (3.8),  Tj^i/2»  ^^+1/2  and  Tl+i/2  are  the  eigenvector 

matrices  of  A,  B  and  C  with  subscript,  e.g.,  in  the  local  ^-direction, 

identifying  that  the  matrix  is  evaluated  at  Qj+i/2,k,l  which  is  a  symmetric  average 

of  Qj,K,L  and  Qj+i,K.L.  expressed  as 

Qj.l,K,L  =  0.5(Qj,k,l-Qj.i,k,l) 
2 

The  m"*  element  of  the  vector  ^j+1/2  is 

Ci  =  'J'(^"i)(gr+gr.i-2a"i) 

where 

g"  =  Smax[o,min(|a'"  i|,Sa,  1)] 
S  =  sign(a"'  1) 

the  difference  of  the  characteristic  variables  aj+y^  is  defined  as 

a^^i  =  r\  (Qj^i,k.l-Qj,k.l) 

2  J*2  ^•^("^J+1,K,L'^Jj,K,l) 


Operator  splitting  and  solution  algorithm 

The  above  TVD  algorithm  can  be  written  in  the  following  operator  form 

LjL^LjAQ"  =  R»  (3.11) 
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where 

=  [I  +  At6jA  +  D,j]" 

=  {l  +  AT[6^C-Re-^6^(r^MJ)]+D,j}'' 
R°  =  -AT[6jE  +  6/  +  6jG-Re-i6jS]"  +  Dg 

Thus,  Eq  (3.11)  can  be  split  into  block  tridiagonal  systems  in  each  coordinate 
direction  and  solved  as  follows 

LjAQ**  =  R"       L^AQ*  =  AQ"       L^AQ"  =  AQ' 

Q"*i  =  Q"  +  AQ" 

in  which  AQ"  and  AQ*  are  intermediate  solutions. 
3.2.2  Numerical  Boundary  Conditions 

Fluid  flow  problems  frequently  involve  a  number  of  fundamentally  different 
types  of  boundaries.  Examples  include  sohd  surface,  far-field,  inflow,  outflow, 
symmetry,  interface  and  singular  boundaries.  In  order  to  obtain  accurate 
solutions,  suitable  boundary  conditions  must  be  imposed  at  each  type  of  boundary. 
The  appropriate  boundary  conditions  for  each  boundary  depend  on  the  nature  of 
the  governing  equations  and  the  domain  transformation.  For  the  thin-layer 
Navier-Stokes  equations,  the  appropriate  boundary  conditions  for  the  sohd 
surface,  inflow,  outflow,  far  field,  and  symmetry  boundaries  were  given  in  Section 
2.6.  In  this  section,  we  discuss  how  these  boundary  conditions  are  implemented 
numerically.  Also,  we  discuss  how  boundary  conditions  are  implemented  at 
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interfaces  where  different  single  grids  of  the  multi-block  grid  meet.  Finally,  we 

discuss  the  averaging  process  imposed  on  the  singular  line. 

Before  discussing  the  details  of  the  boundary  conditions,  we  note  that  for 

the  six-block  topology  used  in  this  study,  all  boundary  conditions  are  treated 

explicitly.  This  method  is  successfully  employed  in  conjunction  with  the  TVD 

implicit  flow  solver,  and  leads  to  a  simple  and  flexible  scheme,  where  boundary 

conditions  can  be  put  into  or  pulled  out  of  a  computer  program  without  disturbing 

the  implicit  algorithm. 

Solid  boundary  conditions 

The  application  of  boundary  conditions  on  the  solid  surfaces  is  simplified 

through  the  coordinate  transformation  since  the  C  =  0  plane  of  each  block 

represent  the  solid  surface  in  the  computational  domain. 

For  viscous  flows  without  injection  or  suction,  the  no-slip  condition  is 

imposed  at  the  sohd  surface  by  setting  the  three  contravariant  velocity 

components  equal  to  zero;  i.e.,  U  =  V  =  W  =  0.  Accordingly,  the  velocity 

components  Ujk,i  =  Vj     =  Wj^j  =  0  for  a  fixed  grid  system. 

Since  viscous  terms  have  a  negligible  effect  on  the  pressure  gradient  at  the 
wall  in  the  direction  normal  to  the  wall,  the  pressure  distribution  along  the  solid  ' 
surface  is  obtained  from  the  three  transformed  momentum  equations  with  the 
viscous  terms  neglected.  By  using  the  metric  identities,  the  viscous  normal 
momentum  equation  for  a  fixed  grid  system  is  given  by  Anderson  et  al.  [26] 
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-  p  W(C  —  +  C  —  +  C  — ) 


in  which 


Since  U  =  V  =  W  =  Oon  the  wall,  the  above  equations  give 

(Vx-^5yCy-^WPt-^(nxCx^nyCy+Ti,Op,-^(CxCx-^CyCy  +  C,gPj  =  0 

This  equation  is  solved  at  the  solid  surface  by  using  second  order  forward 
differences  in  the  C  direction  and  second  order  central  differences  in  the  $  and  ti 
directions  to  obtain  the  surface  pressure. 

The  boundary  condition  imposed  on  the  surface  density  is  dp /da  =  0  at 
the  wall.  The  total  energy  then  is  computed  from  the  known  velocities,  pressure, 
and  density  at  the  soHd  surface  by  using 

Far-field  boundary  conditions 

In  this  study,  the  far-field  boundary  conditions  at  L  =  LMAX  are  set  by 
letting  the  variables  in  the  unknown  vector  Q  equal  to  the  undisturbed  free- 
stream  conditions,  i.e., 
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%^MAX  ~    ^o.         ^J,K,LMAX  ~    ^-  ^J,K,LMAX  ~ 

Pj,K,LMAX  ~     P»  ^J,K,LMAX  ~ 

Inflow/outflow  boundary  conditions 

The  inflow  and  outflow  boundaries  are  special  cases  of  far-field  boundaries 
with  boundary  conditions  that  should  be  specified  according  to  characteristic 
theory.  Here,  the  technique  described  in  the  original  ARC3D  is  used  for  the 
present  study.  At  the  inflow  boundary,  all  variables  are  set  equal  to  their  subsonic 
free-stream  conditions.  At  the  outflow  boundary,  all  variables  are  extrapolated  by 
first-order  linear  extrapolation. 
Interface  boundary  conditions 

An  interface  boimdary  behaves  both  as  an  inflow  boundary  and  as  an 
outflow  boundary.  For  example,  if  flow  enters  a  block  through  an  interface 
boundary,  then  that  interface  behaves  as  an  inflow  boundary  as  far  as  that  block  is 
concerned;  but,  for  the  block  where  the  flow  leaves  it,  that  interface  is  viewed  as 
an  outflow  boundary.  Here,  a  simple  weighted  interpolation  technique  is  used  to 
obtain  the  information  of  unknown  Q  vector  at  the  interface  boundaries  between 
neighboring  blocks.  For  instance,  at  an  interface  boundary  point  of  Blocks  1  and 
2,  the  unknown  variable  q  can  be  estimated  from  the  known  quantities  q^  and  qj 
at  interior  point  of  Blocks  1  and  2  by  the  following  formula 

(Sj  +S2) 
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where  and  are  the  distances  between  the  interior  points  and  the  boundary 
point.  This  formula  assumes  that,  the  gradient  of  q  in  the  local  one-dimension 
direction  s  is  constant;  i.e., 

(\  -  (%      or  -  " 


ds         d&  Sj  s. 


'2 


Singular  boundary  conditions 

As  shown  in  Figure  3.2(a),  the  computational  plane  AEJF  (J  =  1  plane)  of 
Block  1  degenerates  into  a  single  line  in  physical  space.  In  other  words,  each  K 
line  on  the  J  =  1  plane  maps  onto  a  single  point  in  physical  space.  For  each  grid 
point  along  a  specified  K  line  on  this  plane,  the  boundary  conditions  are  obtained 
by  first  using  a  linear  extrapolation  based  on  the  points  at  J  =  2  and  3,  as  shown 
below 

%K,L~  %K,L  ^  %K,L~  q3,K,L 
^12 

where  s^j  and  Sjj  are  the  distances  between  the  points  at  J  =  1  to  J  =  2,  and 
J  =  2  to  J  =  3,  respectively.  Then,  the  average  value  for  all  the  grid  points  which 
along  the  specified  K  line  on  this  J  =  1  plane  will  be  set  as  the  boundary 
conditions  for  each  of  these  grid  points;  i.e.,  the  unknown  variable  q  at  the 
singular  point  is  given  by 
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Symmetry  boundary  conditions 

The  boundary  conditions  at  symmetry  planes  are  treated  by  using  reflected 
image  grid  points  beyond  the  symmetry  boundary.  As  shown  in  Figure  3.4,  the 
z  =  0  surface  in  physical  space  is  also  a  ti  =  0  plane  in  computational  domain, 
which  is  assumed  to  represent  a  symmetry  plane  with  another  plane  added  as  the 
mirror  image  of  the  plane  next  to  the  symmetry  plane.  For  instance,  if  K  =  2  is 
the  symmetry  plane,  then  K  =  1  is  the  added  mirror  image  of  K  =  3  plane,  and 
the  flow  conditions  at  K  =  1  plane  can  be  specified  to  enforce  symmetry  as  shown 
below 

=    Uj3,L  Vj  J  L  =    Vj3  L  Wj  i  L  =  "Wj^l 

Pj,l,L  ^     Pj,3,L  ^J,1,L  ^  ®J,3,L 
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Figure  3.4  Symmetric  boundary  condition 


CHAPTER  IV 
GRID  GENERATION 


In  order  to  use  the  finite-difference  method  presented  in  Chapter  III  to 
obtain  solutions  for  governing  equations  formulated  in  Chapter  11,  it  is  necessary 
to  replace  the  spatial  domain  of  the  problem  being  studied  by  a  finite  number  of 
grid  points.  The  process  of  replacing  a  spatial  domain  by  a  system  of  grid  points 
is  referred  to  as  grid  generation.  In  this  chapter,  we  describe  in  detail  the 
algebraic  grid  generation  technique  used  in  this  study. 

When  a  structured  grid  system  is  used  with  a  finite-difference  method  to 
obtain  solutions  to  fluid  flow  problems,  a  number  of  desirable  conditions  can  be 
listed  to  help  guide  the  grid  generation  procedure,  such  as  the  following:  (1)  For 
computational  efficiency,  the  total  number  of  grid  points  in  the  grid  system  should 
be  kept  to  the  minimum  needed  for  the  finite  difference  method  to  yield  the 
desired  accurate  solutions.  (2)  For  convenience  of  assigning  boundary  conditions, 
one  set  of  grid  lines  always  should  coincide  with  the  boundary  of  the  spatial 
domain  regardless  of  the  geometric  complexity.  (3)  For  ease  of  computing 
derivative  boundary  conditions,  grid  lines  should  intersect  boundaries 
orthogonally.  Orthogonality  is  also  important  near  solid  boundaries  for  turbulent 
flow  problems,  because  the  implementation  of  turbulence  models  often  are  based 
on  information  in  the  normal  direction  [27].  (4)  For  convection  dominated  flows 
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or  when  the  thin-layer  Navier-Stokes  equations  are  used  to  study  such  flows,  one 
set  of  grid  lines  should  be  aligned  with  the  flow  direction.    (5)  The  spacings 
between  grid  points  should  change  slowly  from  a  region  where  grid  points  are 
concentrated  to  a  region  where  grid  points  are  sparsely  distributed. 

For  flow  past  a  simple  geometrical  configuration,  if  we  know  the  locations 
of  important  solution  gradients,  it  is  possible  to  generate  a  good  single  grid.  But 
for  complicated  three-dimensional  flows  within  geometrically  complex  domains, 
the  flow  behavior  may  vary  considerably  from  one  region  to  another.  Therefore, 
it  is  usually  quite  difficult  to  obtain  a  reasonable  single  grid  with  the  entire 
physical  region  transformed  to  a  single  boundary-fitted  coordinate  system  that 
would  satisfy  all  of  the  above  conditions  at  every  part  of  the  spatial  domain.  For 
such  problems,  it  is  often  more  convenient  to  divide  the  entire  spatial  domain 
into  several  blocks  with  each  block  having  a  different  boundary-fitted  coordinate 
system.  With  the  domain  partitioned,  it  is  much  easier  to  generate  a  number  of 
single  grids,  each  of  which  satisfies  the  above  five  conditions.  These  single  grids 
can  then  be  patched  together  by  joining  at  the  interface  boundaries  to  form  the 
overall  grid  system  for  the  entire  flow  domain. 

As  mentioned  in  Section  3.1,  a  multi-block  grid  topology  was  developed  for 
the  wing-fuselage  configuration.  The  flow  domain  was  divided  into  six  blocks  and 
each  block  was  partially  bounded  by  a  solid  surface.  The  challenge  here  is  to 
generate  a  quality  grid  system  for  each  of  the  six  blocks  and,  the  grid  system  in 
the  different  blocks  must  meet  each  other  continuously  at  least  up  to  the  first 
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derivative  at  the  interface  boundaries.  An  overview  of  the  grid  generation 
procedure  used  in  this  study  is  as  follows:  First,  we  use  the  algebraic  surface 
generation  method  developed  by  Shih  et  al.  [16]  to  regenerate  a  smoother, 
contiguous  outer  surface  for  Block  4~the  wing  portion  of  the  domain  that  is  most 
complicated.  Then,  a  one-dimensional  solution-adaptation  algorithm  developed  in 
this  study  is  used  to  cluster  grid  points  for  both  the  inner  and  outer  surfaces  of 
Block  3  and  Block  4  in  regions  where  the  measured  pressure  gradients  on  the 
wing  surface  are  large  and  scattering  them  elsewhere.  Finally,  the  Two-Boundary 
Method  is  used  to  generate  each  single  grid  of  the  multi-block  grid.  These 
methods  are  chosen  because  of  their  high  efficiency  and  ability  to  provide  very 
precise  control  over  the  distribution  of  grid  points  in  the  spatial  domain. 

In  the  following  section,  we  assume  that  the  inner  and  outer  surfaces  of 
each  block  are  known  and  discuss  the  theoretical  background  of  generating 
volume  grids  by  using  the  Two-Boundary  Method.  In  the  section  which  follows, 
we  discuss  the  details  of  the  surface  generation  method  of  Shih  et  al.  and  the  one- 
dimensional  solution-adaptation  algorithm. 

4.1  Volume  Grid  Generation 

In  this  section,  the  theoretical  background  of  the  Two-Boundary  Method  is 
described.  This  discussion  follows  closely  those  of  Yang  and  Shih  [15]  and  Shih  et 
al.  [16]. 
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The  Two-Boundary  Method  is  intended  for  problems  in  which  it  is  only 
necessary  to  map  correctly  two  arbitrary  shaped  opposing  boundaries  of  the 
spatial  domain.  Also,  we  should  note  that  these  two  selected  opposing  boundaries 
of  the  spatial  domain  must  not  intersect  each  other  at  any  point. 

The  first  step  in  generating  volume  grids  by  the  Two-Boundary  Method  is 

to  define  the  coordinate  transformation  between  the  spatial  domain  and  the 

computational  domain,  such  that 

r  =  r(i,r],0 

(4.1) 

=  xU,Ti,C)j+yan,C)k^z(e,Ti,C)i 

In  the  above  equation,  r  represents  the  position  vector  for  the  x-y-z  coordinate 
system  of  the  spatial  domain;  j,  k,  and  1  are  unit  vectors  pointing  in  the  x-,  y-,  and 
z-  directions,  respectively;  i,  x]  and  C  comprise  the  boundary-fitted  coordinate 
system  of  the  computational  domain  which  were  mentioned  in  Section  3.1. 

The  second  step  is  to  specify  the  inner  and  outer  surfaces  which  must  be 
mapped  correctly  as  surfaces  1  and  2,  and  map  them  to  the  coordinate  planes  at 
C  =  0  and  C  =  1  in  the  $-ti-C  coordinate  system,  respectively;  i.e., 

=  r($,Ti,C=0)  =  Ri(5,n) 
'        =  r($,ti,C  =  l)  =  R2a,r\) 

Here,  Rj  is  the  position  vector  of  surface  1  which  has  the  spatial  coordinates  X^,  . 
Yj  and  Zj;  R2  is  the  position  vector  of  surface  2  which  has  the  spatial  coordinates 
X2,  Y2  and  Zj  in  the  x-y-z  coordinate  system.  The  remaining  four  boundary 
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surfaces  3,  4,  5  and  6,  are  mapped  to  coordinate  planes  ti=0,  t|=1,  $=0  and 
5  =  1,  respectively. 

After  specifying  the  two  boundary  surfaces,  the  next  step  is  to  define  curves 
which  connect  these  two  boundaries  by  transfinite  interpolation  techniques.  In 
order  to  implement  the  derivative  boundary  condition  accurately,  the  technique 
based  on  Hermite  transfinite  interpolation  is  used.  This  technique  allows 
specification  of  derivatives  at  the  end  points  of  curves,  ensuring  good 
orthogonality  at  the  boundary  surfaces. 

By  using  the  Hermite  interpolation  to  generate  connecting  curves  between 
surfaces  1  and  2,  the  resulting  cubic  curves  have  the  following  form  as  [16] 

r($,ti,C)  =  Ri(e,n)hi(C)+R2($,Ti)h2(C) 

The  functions  hi(C),  hjCC),  h3(C)  and  h^iC)  are  blending  functions  which 
connect  two  points,  one  on  each  surface  having  the  same  $  and  x]  values.  They 
are  constrained  by  the  following  expressions: 

h,(C=0)  .  1,     h,({.l)  =  0,      ^''l      =  0.  ' =  0 


h,(c=o)=o,  h3(CM)  =  o.  fie^.o 
h,(c.o)  =  0,   h,(c=i)  =  0,  =  o,  =  i 


=  0) 

dh,(C  - 

0) 

dh,(C  = 

0) 

0) 

ac 

ah/c  = 

=1) 

ac 

ah^cc  = 

1) 

ac 

ah3(c  = 

1) 

ac 

ah4(c  = 

1) 

ac 
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With  these  constraints,  hj(C),  h2(C),  and  h4(C)  become 

hi(C)  =  2C3-3C2  +  1 


(4.3) 

h,(0  =  C'-2C'  +  C 


h4(0  =  C'-C^ 

ar(e,Ti,C=0)        ar($,Ti,C  =  l) 
Now,  we  need  to  choose  values  for  rr  and  so  that 

the  connecting  curves  will  intersect  surfaces  1  and  2  orthogonally.  For  surface  1, 
orthogonality  will  occur  when  the  vector  tangent  to  the  connecting  curve  (T^i)  is 
parallel  to  a  unit  vector  normal  to  surface  1  (n^j).  With  an  additional  term 
added  to  control  the  shape  of  the  connecting  curve,  the  derivative  term  is  given  as 

To  =  =  K,(5,n)n,,  =  K,(^,n)  g^g  (4.4) 

'  aT'' air' 

Similarly,  the  derivative  term  for  surface  2  is  defined  as 

To  =  =  K,(J,n)n„  =  K,({,n)  ^  (4.5) 

ail ' 

The  coefficients  Ki($,ti)  and  K2($,n)  appearing  in  the  above  equations  are  known 
as  the  K-factors  which  may  be  constants  or  variables  and  are  chosen  by  numerical 
experiments  so  that  no  overlapping  of  the  connecting  curves  takes  place  in  the 
interior  of  the  spatial  domain. 
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The  desired  cubic  connecting  curves  based  on  the  Hermite  interpolation 
can  be  obtained  by  substituting  Eqs.  (4.3)  to  (4.5)  into  Eq.  (4.2). 

Having  mapped  the  x-y-z  coordinate  system  onto  the  l-r\-C  coordinate 
system,  we  now  proceed  to  discretize  the  ^-r]-(  coordinate  system  by  some  equally 
spaced  grid  points.  The  locations  of  these  grid  points  are  given  by  the  ordered 
triples  ($j,  tiK,  ^l)  where 


5j 

=  (J-1)A$ 

J  = 

1,2,-,JMAX 

=  (K-i)An 

K  = 

1,2,-,KMAX 

(4.6) 

=  (L-l)AC 

L  = 

1,2,  -LMAX 

and 

=   1  ,       An  =   1  ,       AC  =   -  

JMAX-1  KMAX-1  LMAX-1 

After  the  locations  of  grid  points  in  the  computational  domain  are  known 
we  can  obtain  the  locations  of  the  grid  points  in  the  x-y-z  coordinate  system  by 
substituting  Eqs.  (4.3)  to  (4.6)  into  Eq.  (4.2). 

At  this  point,  we  need  to  consider  the  physics  of  our  problem  for  which  the 
volume  grid  points  are  generated.  In  order  to  get  accurate  solutions  for  the 
viscous  transonic  flow  problems,  we  need  to  resolve  the  viscous  sublayer.  Hence, 
the  grid  points  should  be  highly  concentrated  in  the  direction  normal  to  the  solid 
surface.  Such  a  grid  spacing  control  can  be  achieved  by  the  use  of  stretching 
functions,  which  is  an  efficient  way  to  redistribute  the  grid  points. 
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Several  different  stretching  functions  have  been  investigated  in  details  [5, 
39,  40],  and,  in  comparison  based  on  the  truncation  error  analysis  [40],  a 
hyperbolic  tangent  stretching  has  the  better  overall  distribution  inside  and  outside 
the  viscous  sublayer.  Therefore,  this  stretching  function  is  used  to  refine  the  grid 
spacings  near  the  solid  boundary. 

The  hyperbolic  tangent  stretching,  e.g.  [5],  is  constructed  by  letting  the 
normalized  arc  length  ST  vary  from  0  to  1  as  L  varies  from  1  to  LMAX,  where 
LMAX  is  the  total  number  of  grid  points  along  the  given  arc  length.  Then, 
specify  the  normalized  values  of  the  minimum  grid  spacing  next  to  the  soUd 
boundary  as  DSO  and  the  ratio  of  last  grid  spacing  DSM  at  each  end  of  arc,  and 
define 


A  = 


DSM 


^J  DSO 

B  =  (LMAX-lVDSMxDSO 
The  next  step  is  to  find  the  distribution  factor  W  by  using  the  Newton- 

Raphson  method  to  solve  the  following  nonlinear  equation 

sinhW  1 


W  B 

Then,  the  normalized  hyperbolic  tangent  arc  length  distribution  is  given  as 

Si  =  0,  Slmax  =  1.  and    S^  =        "\  L  =  2,3,-,  LMAX -1 

A  +  (1-A)Ul 

where 

tanh[W( — tlL_-i)] 

u,  =  Uu  LMAX-1  2^\ 

2  W 
tanh(^) 
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The  distribution  of  grid  points  in  the  x-y-z  coordinate  system  along  each  arc 
can  be  obtained  by  rescahng  the  normaUzed  arc  length  distribution  back  to  the 
original  scale  and  interpolating  by  a  simple  and  computationally  efficient  Aitken's 
algorithm  [41], 

Finally,  the  last  step  is  to  smooth  discontinuities  in  the  grid  lines  and 
reduce  any  rapid  spacing  changes  between  neighboring  grid  points  in  the  spatial 
domain.  The  method  of  smoothing  used  here  is  to  replace  the  coordinates  at  a 
interior  grid  point  by  a  weighted  average  of  itself  and  the  two  adjacent  points 
along  a  grid  line,  i.e. 

''j,K,L  "  '^J,K,L'^"(^J+1,K,L~2Xjk,l*^J-1,K,l) 

yj,K,L  =  yj.K,L+"(yj.i,K,L-2yj,K,L-^yj-i,K.L)  (4.7) 

in  which  the  value  of  w  used  is  0.25  to  keep  the  smoothing  procedure  stable. 

4.2  Surface  Generation 

An  essential  aspect  of  generating  grid  networks  for  the  six-block  wing- 
fuselage  flow  problem  by  using  the  Two-Boundary  Method  is  that  the  two  suitable 
boundary  surfaces  of  each  block  must  be  constructed  first.  The  two  boundary 
surfaces  of  the  six-block  grids  are  first  obtained  from  a  previous  study  [25].  The 
outer  surface  of  Block  4  is  regenerated  by  using  the  three-dimensional  Bi- 
Directional  Hermite  Interpolation  Method  [16],  while  the  four  boundary  curves  of 
these  two  outer  surfaces  are  modified  from  the  previous  discontinuous  surfaces  to 
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ensure  the  obtaining  of  contiguous  grid.  Afterwards,  grid  points  on  both  the  inner 
and  outer  surfaces  of  Block  3  and  Block  4  are  redistributed  by  a  one-dimensional 
solution-adaptive  method.  Finally,  grid  points  on  both  the  iimer  and  outer 
surfaces  of  Block  2  and  Block  5  are  redistributed  according  to  the  distribution  of 
grid  points  on  Block  3  and  Block  4,  respectively,  by  using  Aitken's  algorithm  to 
ensure  contiguous  and  smooth  surfaces. 
4.2.1  Three-Dimensional  Bi-Directional  Hermite  Interpolation 

In  the  following  discussion,  we  will  describe  the  theoretical  background  of 
generating  the  surface  which  is  mapped  onto  the  coordinate  plane  located  at 
C  =  1  in  the  computational  domain  by  the  Bi-Directional  Hermite  Interpolation 
Method. 

To  generate  a  three-dimensional  surface  by  using  the  Bi-Directional 
Hermite  Interpolation,  the  first  step  is  to  seek  a  coordinate  transformation  and 
represent  the  position  vector  r  as  expressed  in  Eq.  (4.1). 

For  surface  generation,  all  four  boundary  curves  need  to  be  mapped 
correctly  to  the  coordinate  lines  in  the         coordinate  system.  Hence,  the 
second  step  is  to  specify  curves  1,  2,  3,  and  4  as  coordinate  lines  n  =  0,  ti  =  1, 
5=0  and  1  =  1,  respectively;  then  represent  the  position  vector  of  the  boundary 
curves  in  the  x-y-z  coordinate  system  as  Rj,  R2,  R3,  and  R4,  i.e., 

Ri  =  ra,Ti=0,C  =  l)  =  Rj(0       R2  =  r($,Ti=l,C  =  l)  =  R,a) 

R3  =  r($=0,n,C  =  l)  =  RjCn)       R4  =  r($=l,Ti,C  =  l)  =  R4(ti) 
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where  the  coordinates  of  these  curves  in  the  x-y-z  coordinate  system  are  expressed 
as  Xj,  Yj,  Zj,  and  i  =  1,  2,  3,  4. 

Once  the  four  boundary  curves  have  been  specified,  the  next  step  is  to 
correctly  map  the  boundary  curves  1  and  2  (which  must  not  touch  each  other  at 
any  point)  and  define  curves  connecting  these  two  boundary  curves  by  the 
following  Hermite  transfinite  interpolation  expression: 
rc(U,C  =  l)  =  R,(£)h,(n)+Rj(?)hj(ti) 

.M!:aiMli)h3(n)*M!iaiMzi)h,(n)  ^''^^ 

With  the  blending  functions  hi(ti),  hjCn),  h3(Ti)  and  h4(n)  which  connect  two 
points,  one  on  each  curve  having  the  same  $  and  C  values  expressed  as 

h,(n)  =  2n3-3Ti2  +  i 

h,(r))  =  -2n3^3ri' 

hgCn)  =  Ti^-2Ti'+n 
Vn)  =  n'-n' 

When  the  above  step  is  completed,  similar  to  the  Two-Boundary  Method 
described  in  Section  4.1,  we  correctly  map  curves  1  and  2  in  the  x-y-z  coordinate 
system  to  coordinate  lines  x]  =  0  and  n  =  1  in  the  i-r]-C  coordinate  system. 
However,  curves  3  and  4  may  in  general  have  not  been  mapped  correctly  to 
coordinate  lines  $  =  0  and  $  =  1.  Instead,  straight  lines  3'  and  4'  have  been 
generated  in  the  spatial  domain,  and  it  is  these  curves  which  have  been  mapped 
to  the  coordinate  lines  $  =  0  and  $  =  1  in  the  computational  domain. 
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In  order  to  correctly  map  curves  3  and  4  in  the  x-y-z  coordinate  system  to 
coordinate  lines  i  =  0  and  $  =  1  in  the  i-r\-C  coordinate  system,  we  need  to 
adjust  curves  3'  and  4'  by  defining  Ar  such  that 

ra,ii,C  =  l)  =  r;2($,Ti,C  =  l)  +  Ara,n,C  =  l)  (4.9) 

will  map  curves  3  and  4  to  $  =0  and  5  =  1,  respectively. 

Here,  Ar  is  given  by  the  following  Hermite  interpolation  expression: 
Ar($,n,C  =  l)  =  [r(5=0,n,C  =  l)-ri5($=0,n,C  =  l)]h5(O 

+  [r($  =  i,n,C  =  1) -ri^a  =  i,n,C  =  l)]h,(0 

,  ,ar(g=o.n,C  =  i)  gr;,($=o.n.C  =  i)^^^^^  (4-io) 


^  rar($  =  i,Ti,c  =  i)_^i2a=i>n>C  =  i)^^ 
^      8$       "       al  ^^^^^ 


where 


ar(,(g  =  i,n,C  =  i)  ^  ar($  =  i,n=o,c  =  i).  ,,,,ar($  =  i,Ti  =  u  =  i)>, ,  , 
a$  al        ^^^^   di  


(4.11) 


(4.12) 
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and 


(4.13) 


e-e 


We  still  need  to  evaluate  the  partial  derivative  terms  appearing  on  the 
right  hand  sides  of  Eqs.  (4.11)  and  (4.12).  Since  the  only  information  we  have 
about  the  surface  is  the  four  boundary  curves,  these  boundary  curves  are  used  to 
derive  expressions  for  the  first-order  derivative  terms. 

We  shall  use  the  method  described  by  Shih  et  al.  [16]  to  evaluate  the  first- 
order  derivative  terms.  For  instance,  to  derive  these  terms  along  curve  1  on 
which  n  =  0,  C  =  1,  and  $  varies  between  0  and  1,  this  method  involves  the 
following  major  steps: 

(1)      Determine  normal  vector  at  each  end  of  curve  1. 

The  vector  perpendicular  to  the  plane  which  formed  by  the  two  vectors 
tangent  to  curves  1  and  3  at  5  =  0  and  n  =  0  is  given  by: 


Similarly,  the  vector  perpendicular  to  the  plane  which  formed  by  the  two  vectors 
tangent  to  curves  1  and  4  at  $  =  1  and  n  =  0  is  given  by: 


N, 


aR,($=o)  aR3(ti=o) 


X 


aR^($  =  i)^aR,(n=o) 

dr\ 
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(2)  Linearly  interpolate  between  two  normal  vectors. 

In  this  step,  we  linearly  interpolate  between  the  two  normal  vectors 
obtained  in  the  above  step  and  get  the  following  vector  function 

N^i  =  (l-ONi3^5N„ 

(3)  Determine  a  vector  function  tangent  to  the  surface. 

Since  there  are  an  infinite  number  of  surfaces  that  can  pass  through  any 
given  curve,  a  strategy  adopted  here  is  to  select  the  surface  which  will  pass 
through  curves  1,  3,  and  4  as  the  approximating  surface  that  partial  derivatives 
with  respect  to  ti  need  to  be  evaluated  on. 

Here,  let  T^j  be  a  vector  function  tangent  to  curve  1  and  given  by  the 
expression  as 

"  ^  ~dr 

Thus,  the  unit  vector  function  tangent  to  the  surface  at  curve  1  can  be 
approximated  by 

(4)  Determine  expressions  for  the  first  order  derivatives.  , 

Since  we  desire  for  grid  lines  on  the  surface  to  be  perpendicular  to  the 
boundary  curve  1  for  this  case,  the  vector  function  tangent  to  the  n -direction,  T^^ 
must  be  parallel  to  the  unit  vector  function,  e^^,  derived  in  step  3.  In  a  manner 
similar  to  the  discussed  derivation  of  the  Two-Boundary  Method,  the  desired 
expressions  for  the  first  order  partial  derivatives  can  be  given  by 

\t  -  —  
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Here,  the  Kj-factor  is  also  chosen  by  numerical  experiments  to  ensure  the 
creation  of  a  smooth  surface  and  to  prevent  overlapping  of  the  connecting  curves 
in  the  interior  of  the  surface.  The  second  order  derivative  terms  in  Eqs.  (4.11) 
and  (4.12)  are  set  equal  to  zero. 

The  derivation  of  partial  derivative  terms  along  the  other  three  boundary 
curves  is  similar  to  the  procedure  described  above.  Thus,  the  desired  surface  in 
the  x-y-z  coordinate  system  can  be  generated  by  substituting  Eqs.  (4.10)  to  (4.13) 
into  Eq.  (4.9). 

4.2.2  One-Dimensional  Solution-Adaptation 

As  mentioned  earlier,  for  computational  efficiency  and  accuracy,  the  grid 
points  should  be  concentrated  in  regions  of  sharp  gradients  or  small  length  scale 
phenomena  and  sparsely  distributed  elsewhere.  In  this  investigation,  the  viscous 
sublayer  is  resolved  by  the  stretching  functions  as  described  in  Section  4.1.  The 
high  pressure  gradients  along  the  streamwise  direction  on  the  wing  surface  is 
resolved  by  a  solution-adaptive  scheme,  by  which  the  physics  of  the  flow  problem 
dictate  the  distribution  of  the  grid  points.  This  algorithm  is  adopted  in  this  study 
because  it  has  been  shown  to  provide  good  grid  resolution  for  complex  transonic 
flow  fields  [42-44]. 

To  develop  a  solution-adaptive  grid  generation  scheme  using  variational 
principle,  a  functional  in  each  coordinate  direction  is  defined  and  minimized  to 
result  in  a  set  of  partial  differential  equations  which  govern  the  grid  point  spacing. 
In  many  instances,  adaptation  is  required  in  only  one  coordinate  direction.  For 
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this  case,  the  partial  differential  equation  reduces  to  an  ordinary  differential 
equation.  Because  of  the  simplicity  of  the  one-dimensional  variational  formula,  it 
is  employed  in  the  present  study  and  described  in  the  following  sequence.  First,  a 
functional  to  be  minimized  by  variation  principle  and  a  weighting  function 
associated  with  the  distribution  of  grid  points  are  defined.  Afterwards,  the 
algorithm  of  using  a  defined  weighting  function  to  create  a  solution-adaptive  grid 
generation  scheme  is  described. 
Variational  principle  and  functional 

In  the  following  discussion,  the  grid  points  along  a  Une  are  specified  by  the 
successive  integer  values  of  $(S),  where  $(S)  =  J;  the  maximum  value  of  $(S)  is 
equal  to  the  total  number  of  grid  points  on  the  line,  i.e.,  JMAX,  and  A$(S)  =  1. 
The  nonuniform  grid  points  distribution  in  the  spatial  domain  can  be  considered 
as  a  transformation  from  a  uniform  distribution  in  $(S)-space  to  the  arc  length 
S($)  along  a  given  curve. 

In  order  to  formulate  a  one-dimensional  solution-adaptive  grid  generation 
scheme,  a  functional  to  be  minimized  must  be  defined  first.  Here,  this  functional 
is  defined  as  the  least-squares  of  the  desired  grid  property  in  the  ^  coordinate, 
G($),  over  a  positive  weighting  function  W($)  which  is  associated  with  the 
distribution  of  grid  points;  and  is  given  by 


In  the  above  functional,  when  the  weight  function  W($)  is  small,  the  grid  property 
G($)  must  also  be  small  in  order  to  minimize  the  integral. 


(4.14) 


63 

Since  we  desire  a  small  weighting  function  to  correspond  to  a  large  grid 
spacing,  the  grid  property  G($)  is  chosen  to  be  the  grid  point  density  and 
expressed  as 

G(0  =      =  ^  =  ^  (4.15) 

where  Sj  represents  the  grid  spacing,  and  a  small  value  of     corresponds  to  a 
large  grid  density.  By  using  Eq.  (4.15),  the  functional  expressed  by  Eq.  (4.14)  can 
now  be  written  as 

I  =  [[  J!-fdS  (4.16) 

The  problem  of  minimizing  the  above  integral  can  be  solved  by  a  variation 
principle  and  the  solution  is  the  Euler-Lagrange  equation  associated  with  that 
integral.  The  Euler-Lagrange  equation  is  given  by 

3$  ds 

where 

F  =  (4.18) 
W(0 

Substituting  Eq.  (4.18)  into  Eq.  (4.17)  and  simplifying  yields  the  following  ordinary 
differential  equation 

ds  W($) 

which  can  be  integrated  to  obtain  the  governing  equation  as 

W($)Sj  =  constant  (4.20) 
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The  discrete  form  of  the  above  governing  equation  is  given  by 

WjASj  =  constant  (4.21) 

where 

ASj  =  Sj,i-Sj  J  =  1,2,-,JMAX-1 

Weighting  function 

To  make  the  method  solution-adaptive,  the  weighting  function  appearing  in 
Eq.  (4.21)  should  be  related  to  the  solution  of  the  flow  problem  and  have  the 
effect  to  control  the  grid  spacing.  For  this  flow  problem,  it  is  chosen  as  a  function 
of  the  pressure  gradient.  However,  letting  the  weighting  function  be  equal  only  to 
the  pressure  gradient  causes  two  problems.  First,  when  the  gradient  is  zero,  Eq. 
(4.21)  is  singular.  Second,  when  the  gradient  is  very  large,  almost  all  of  the  grid 
points  will  be  clustered  in  the  vicinity  of  the  large  gradient  so  that  the  grid  system 
will  not  be  smooth.  Hence,  a  proper  weighting  function  must  have  terms 
controlling  both  grid  smoothness  and  grid  adaptation.  In  similar  to  a  previous 
study  [44],  the  weighting  function  used  here  has  the  form  of 

Wj  =  1+aff  (4.22) 
which  does  control  the  grid  smoothness  and  adaptation.  With  a  limited  number 
of  grid  points  available  for  use,  a  proper  choice  of  the  unknown  constants  a  and  p 
can  result  in  a  good  grid  adaptive  to  the  normalized  magnitude  of  pressure 
gradient  f,. 

The  nonnegative  constants  a  and  p  shown  in  Eq.  (4.22)  can  be  determined 
by  specifying  the  minimum  and  maximum  grid  spacings,  AS^^^  and  AS„i„ ,  along  a 
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coordinate  line  for  (fj)„ax  =  1-0  and  (fj)^^  =  0.0,  respectively.  Accordingly,  Eqs. 
(4.21)  and  (4.22)  give 

AS  .  / 

a  =  "^-1 
AS 

mm 

Note  that  proper  choices  of  AS^i„  and  AS^^x  depend  on  the  physics  of  a  problem 
as  well  as  on  the  number  of  grid  points  available.  Furthermore,  we  have  from 
Eqs.  (4.21)  and  (4.22)  the  relation 


JMAX-1 


=  constant 


J-i  1+afJ 

in  which  Sj  is  the  total  length  of  a  given  curve.  Since  Wj  is  a  maximum  when  AS 
is  a  minimum,  AS„i„  can  be  expressed  as 

JMAX-1  I  (4-23) 


J 


=1  Uaf; 


P 


Accordingly,  the  remaining  unknown  constant  p  can  be  determined  from  the 
above  nonlinear  equation.  Let  ' ;    i  '      '  - 

JMAX-1       ^  c 

with  '  f-'  I. ,  * 

,   ,a\  JMAX-1  -,fP 

g^(P)  =  ^=AS„,  E  [-^ln(f.)] 
°P  J=i  (l+afjV 

then  p  can  be  found  by  the  Newton-Raphson  method  for  g(p)  =  0. 


66 

After  determining  the  two  constants  and  weighting  function,  the  grid  points 
can  be  redistributed  according  to  Eq.  (4.21)  to  make  the  product  of  weighting 
function  and  grid  spacing  constant,  and  ensuring  that  the  grid  points  will  be 
distributed  so  that  the  grid  spacing  is  small  where  the  weighting  function  is  large, 
and  vice  versa. 


CHAPTER  V 
RESULTS  AND  DISCUSSION 


In  this  chapter,  we  demonstrate  the  usefulness  of  the  multi-block 
computational  method  being  developed.  This  method  is  indeed  rather  simple  and 
straightforward  in  application;  however,  special  features  of  the  method  have 
resulted  in  excessive  distortions  on  domain  transformation  for  complex 
configuration  flow  problem  and  consequently  make  it  difficult  to  generate  quality 
block  grids  for  accurate  flow  simulations.  In  order  to  gain  further  insight  and 
understanding  on  the  developed  method,  a  number  of  numerical  experiments  with 
emphasis  on  the  effect  of  block  griddings  to  the  solution  accuracy  have  been 
conducted  in  this  study.  As  described  in  Section  2.1,  the  flow  problem  being 
simulated  is  a  transonic  turbulent  flow  of  Mach  0.8  past  a  wing-fuselage 
configuration,  shown  in  Figure  2.1,  at  -3  degree  angle  of  attack.  In  the  solution 
method,  the  flow  field  of  interest  is  properly  divided  into  six  contiguous  blocks, 
Block  l~nose  block.  Blocks  2  and  5~side  blocks.  Blocks  3  and  4~wing  blocks,  and 
Block  6~tail  block. 

As  reported  in  a  previous  study  [25],  the  imposed  grid  topology  and 
boundary-fitted  domain  transformation  did  make  it  very  difficult  to  generate  a 
smooth  and  nonsingular  block  grid  for  the  wing  blocks,  since  four  of  the  six 
bounding  planes  in  the  computational  domain  form  a  nearly  flat  surface  in  the 

67 


68 

spatial  domain.  In  fact,  the  wing  block  grids  generated  and  used  in  flow 
simulation  had  negative  Jacobians  of  transformation  at  many  grid  points  when 
they  were  evaluated  by  finite-difference  approximations.  However,  with  the 
inverse  of  Jacobian  approximated  by  the  grid  cell  volume,  numerical  simulation  of 
the  flow  problem  at  zero  degree  angle  of  attack  had  resulted  in  reasonably 
accurate  prediction  of  shock  wave  location  and  acceptable  agreement  between  the 
measured  and  computed  wing  surface  pressure  distribution.  More  accurate 
solutions  might  have  resulted  if  better  wing  block  grids  had  been  generated  and 
used  in  the  flow  simulation.  For  such  a  complex  configuration  flow  problem,  a 
better  block  grid  can  be  measured  with  the  existence  of  the  number  of  negative 
Jacobians  evaluated  by  the  finite-difference  approximation. 

In  this  study,  the  six  block  grids  are  generated  with  the  grid  generation 
techniques  described  in  the  previous  chapter.  For  adaptive  wing  surfaces  having 
40  grid  points  in  the  chordwise  direction,  the  distribution  of  the  40  points  is 
determined  by  the  one-dimensional  adaptive  gridding  in  which  the  adaptive 
function  ^  in  Eq.  (4.22)  is  provided  with  the  measured  and  interpolated  surface 
pressure  gradient.  The  maximum  and  minimum  grid  spacings  chosen  for  the 
adaptation  are  AS^^x  =  l-5h,  and  AS„i„  =  0.4h;  here  h  represents  the  averaged 
grid  spacing  along  a  grid  line.  In  the  volume  grid  generation  procedure,  grid 
orthogonality  is  not  enforced  for  the  outer  surfaces.  Therefore,  the  Kj-factors  in 
Eq.  (4.5)  are  set  equal  to  zero.  However,  the  K^-factors  in  Eq.  (4.4)  for  the  inner 
surfaces  are  chosen  by  numerical  experiments  to  ensure  nonoverlapping  of  grid 
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lines.  For  each  block,  the  first  grid  spacing  away  from  the  solid  surface  is  set  to 
about  2xlO"*D,  where  D  is  the  largest  cross-section  of  fuselage  (see  Section  2.1), 
which  is  considered  fine  enough  to  resolve  the  viscous  boundary  layer  for 
turbulent  flow  simulations,  and  the  last  grid  spacing  next  to  the  outer  boundary  is 
set  to  ten  times  that  of  the  average  spacing  along  a  given  grid  line. 

Five  different  sets  of  block  grids  have  been  generated  and  provided  to  the 
TVD  thin-layer  Navier-Stokes  code  for  the  flow  simulation.  The  computation  was 
performed  on  a  Cray  Y-MP  system  at  Florida  State  University.  For  the  flow 
problem  studied,  measured  surface  pressure  coefficients  (Cp)  at  eight  wing  cross- 
sections  are  available  for  comparison  with  the  computed  results.  As  shown  in 
Figure  5.1,  these  cross-sections  are  at  the  normalized  distance  0.150,  0.225,  0.305, 
0.405,  0.500,  0.700,  0.850,  and  0.950  measured  from  the  symmetry  plane  of  the 
configuration.  The  measured  Cp  distribution  represented  by  the  symbols  °  and  o 
for  the  lower  and  upper  wing  surfaces,  respectively,  are  shown  in  Figure  5.2  for 
two  representative  cross-sections.  Results  obtained  in  this  study  are  presented 
and  discussed  in  the  following  order.  First  the  importance  of  adaptive  wing 
surface  grids  is  identified  in  Section  5.1.  Afterwards,  the  effect  of  grid  resolution 
to  the  solution  accuracy  is  discussed  in  Section  5.2.  Finally  the  complex  flow 
physics  simulated  are  presented  in  Section  5.3. 
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Uing  tip 


Figure  5.1  Normalized  distances  on  wing  surface 


(a)  d  =  0.305  (b)  d  =  0.700 

Figure  5.2  Measured  Cp  distribution  plot 
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5.1  Effect  of  Solution-Adaptation 

The  first  set  of  six  block  grids,  Grid  set  1,  considered  has  the  dimension  of 
28x25x60,  40x13x60,  40x35x60,  40x35x60,  40x13x60,  and  26x25x60,  respectively,  for 
Block  1  to  Block  6.  The  inner  and  outer  surfaces  employed  for  block  grid 
generations  are  obtained  from  the  previous  study  [25],  except  that  the  outer 
surface  of  Block  4  has  been  modified  to  ensure  the  continuity  of  the  block 
surfaces.  The  method  used  for  the  modification  is  the  three-dimensional  Bi- 
Directional  Hermite  Interpolation  Method  discussed  in  Section  4.2;  the 
intermediate  surface  generated  from  Eq.  (4.8)  is  presented  in  Figure  5.3(a),  and 
the  resulting  surface  obtained  from  Eq.  (4.9)  is  shown  in  Figure  5.3(b).  For  this 
set  of  grid  network,  the  distribution  of  grid  points  on  each  surface  is  not  adaptive 
to  any  solution  field.  The  representative  grid  characteristics  with  grid  points 
concentrated  in  the  leading  and  trailing  edge  area  can  be  observed  from  Figure 
5.4  and  Figure  5.5  for  the  inner  surfaces  of  Block  3  and  4,  respectively. 

The  volume  grid  of  each  block  is  generated  by  using  the  Two-Boundary 
Method  discussed  in  Section  4.1.  Numerical  experiments  have  been  conducted 
and  results  obtained  indicate  that  reasonably  good  grids  can  be  generated  for  each 
of  the  six  blocks  if  the  values  of  K^-factors  defined  in  Eq.  (4.4)  are  chosen  as  0., 
0.,  0.001,  0.0018,  0.,  and  0.,  respectively,  for  Block  1  to  Block  6.  Figure  5.6 
presents  the  generated  volume  grid  of  Block  3  with  cut-offs  for  reference.  An 
assessment  of  the  six  block  grids  generated  shows  that  the  Jacobian  approximated 


(b)  Resulting  surface 
Figure  5.3  Modified  outer  surface  of  Block  4 
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(b)  K  =  1  to  20 
Figure  5.4  Inner  surface  grid  points  distribution  of  Block  3  (Grid  1) 
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(b)  K  =  1  to  20 
Figure  5.5  Inner  surface  grid  points  distribution  of  Block  4  (Grid  1) 
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Figure  5.6  Block  3  -  volume  grids  with  cut-offs  (Grid  1) 
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by  finite-difference  approximations  has  resulted  in  negative  value  at  many  grid 
points  in  the  wing  blocks.  Accordingly,  in  the  flow  simulation  process,  the 
Jacobian  has  been  represented  by  the  grid  cell  volume. 

With  the  maximum  nondimensionalized  time  step  (Axn,a^)  shown  in  Eq. 
(3.11)  set  to  1.0,  the  calculated  Cp  distribution  obtained  in  4000  time  steps  (the 
total  nondimensionalized  time  Xj^,,  =  2750)  for  four  representative  cross-sections 
on  the  wing  surfaces  is  presented  in  Figure  5.7.  From  this  figure  with  the 
computed  Cp  distribution  on  the  lower  and  upper  wing  surfaces  represented  by 
the  dashed  and  solid  lines,  respectively,  we  can  observe  that  although  the 
calculated  data  agree  qualitatively  with  the  experimental  data,  large  discrepancies 
can  be  observed  in  the  leading  edge  area,  around  the  shock  wave,  and  high 
pressure  gradient  regions. 

In  order  to  improve  the  grid  distribution,  the  one-dimensional  solution- 
adaptation  algorithm  described  in  Section  4.2  is  used  to  redistribute  the  grid 
points  on  the  inner  and  outer  surfaces  of  Block  3  and  Block  4.  The  resulting 
adaptive  wing  surfaces  are  shown  in  Figure  5.8  and  Figure  5.9  for  comparison  with 
those  shown  in  Figure  5.4  and  Figure  5.5.  Afterwards,  grid  points  on  the  inner 
and  outer  surfaces  of  Block  2  and  Block  5  are  also  redistributed  according  to  the 
locations  of  grid  points  on  those  two  boundary  surfaces  of  Block  3  and  4  by  the 
Aitken's  algorithm  to  ensure  contiguous  interfaces  and  smooth  surfaces.  Finally, 
the  second  set  of  volume  grids  are  generated.  Numerical  experiments  conducted 
have  resulted  in  the  choice  of  appropriate  values  for  the  Kj-factors  of  0.,  0., 
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Figure  5.7  Cp  distribution  plot  (Grid  1) 


78 


(a)  K  =  20  to  35 


(b)  K  =  1  to  20 
Figure  5.8  Inner  surface  grid  points  distribution  of  Block  3  (Grid  2) 
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0.  0005,  0.0005,  0.,  and  0.,  respectively,  for  Block  1  to  Block  6.  Similar  to  Grid  set 

1,  the  block  grids  generated  have  negative  Jacobians  at  many  grid  points,  which 
could  be  attributed  to  large  truncation  error.  Again,  the  Jacobian  is  approximated 
by  grid  cell  volume  in  the  flow  solver. 

With  AT^a^  set  to  1.0,  the  computed  Cp  distribution  obtained  in  4000  time 
steps  (t^ot  =  2050)  is  shown  in  Figure  5.10.  When  compared  with  Figure  5.7,  one 
can  observe  the  effect  of  adaptive  gridding  which  has  improved  the  agreement  in 
the  high  pressure  gradient  regions.  However,  the  agreement  with  measured  data 
has  deteriorated  in  other  regions,  apparently  due  to  the  use  of  insufficient  grid 
points. 

For  the  two  cases  investigated,  the  Cray  Y-MP  CPU  time  required  to 
march  the  solution  for  1000  time  steps  calculation  is  about  2.95  hovu-s.  The  , 
unsatisfactory  quantitative  agreement  could  be  a  result  of  insufficient  grid 
resolution  on  the  wing  surfaces.  However,  the  overall  agreement  between  the 
measured  data  and  computed  results  is  quite  encouraging. 

5.2  Effect  of  Grid  Refinement 

Three  different  sets  of  refined  block  grids,  Grid  sets  3,  4,  and  5,  have  been 
considered  for  investigating  the  effect  of  grid  resolutions  to  the  accuracy  of  the 
flow  simulation  as  well  as  to  the  efficiency  of  the  solution  method.  The 
dimension  of  the  three  grid  sets  generated  is  summarized  in  Table  5.1.  Also,  the 
maximum  allowable  time  step  At^^^  in  the  flow  simulation,  total 
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(c)  d  =  0.700  (d)  d  =  0.850 


Figure  5.10  Cp  distribution  plot  (Grid  2) 
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Table  5.1  Grid  dimension  and  computational  information  (Grid  set  3  to  5) 


frrid  set 

Grid  set  4 

Rlnrk  1 

5Sx2SxfiO 

55x47x60 

Rlnrk  ? 

1 S7Y?4Y6n 

1 J  /aHUaOU 

Block  4 

79x35x60 

79x69x60 

157x40x60 

Block  5 

79x13x60 

79x13x60 

157x24x60 

Block  6 

51x25x60 

51x25x60 

51x47x60 

^^max 

0.20 

0.20 

0.05 

"^Tot 

1258 

1368 

672 

CPU  hours 

5.72 

8.89 

13.5 

nondimensionalized  time  t^o,,  and  the  Cray  Y-MP  CPU  time  required  per  1000 
time  steps  of  calculation  are  given  in  the  table  for  reference. 

Grid  set  3  is  designed  to  investigate  the  effect  of  grid  resolution  in  the 
streamwise  direction.  The  inner  and  outer  surfaces  for  use  in  the  block  grid 
generation  are  obtained  from  those  of  Grid  set  2  by  halving  the  grid  spacings  in 
$ -direction.  Hence,  the  refined  vdng  surfaces  are  also  adaptive  to  the  measured 
pressure  gradient.  The  volume  grids  of  each  block  are  then  generated  by 
choosing  the  values  of  K^-factors  as  0.,  0.,  0.001,  0.005,  0.,  and  0.,  respectively,  for 
Block  1  to  Block  6.  An  assessment  of  the  refined  block  grids  generated  shows 
that  the  Jacobian  evaluated  by  finite  difference  approximations  is  positive  at  every 
grid  point  by  one  of  the  following  two  smoothing  procedures.  The  first  smoothing 
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procedure  is  described  in  Section  4.1.  The  second  smoothing  procedure  is  given 
by 

X2  =  Xx^  +  (1-X)x^ 

Yi  =  ^yi+(i-^)y3 

Z2  =  A.Zi+(l-X)Z3 

where  A.  is  a  constant  determined  by  numerical  experiments,  and  (xj,  y2,  Z2)  is  a 
point  on  the  grid  Une  that  overlaps  either  the  grid  line  that  passes  through  point 
(xj,  yi,  Zi)  or  the  point  (X3,  yj,  Z3).  The  Cp  distribution  calculated  at  the  4000th 
time  step  (tt^,  =  1258)  and  the  experimental  data  at  four  different  cross-sections 
are  presented  in  Figure  5.11.  Observe  that  the  predicted  Cp  distribution  and 
shock  wave  location  on  the  lower  wing  surface  are  in  much  better  agreement  with 
measured  data  than  those  shown  in  Figure  5.10;  however,  the  predicted  result  on 
the  upper  wing  surface  has  little  improvement.  Note  that  the  maximum  time  step 
allowed  by  the  flow  solver  has  been  reduced  from  1.0  to  about  0.20. 

In  order  to  investigate  the  effect  of  grid  resolution  on  the  wing  spanwise 
direction.  Grid  set  4  is  generated  from  Grid  set  3  by  halving  the  grid  spacing  in 
the  wing  spanwise  direction.  As  expected  the  Jacobian  evaluated  by  finite 
difference  approximations  is  positive  at  every  grid  point.  The  Cp  distribution 
obtained  after  4000  time  steps  of  computation  (x^^^  =  1368)  are  shown  in  Figure 
5.12,  for  comparison  with  those  of  Figure  5.11.  It  is  clear  that  the  refined  grid 
resolution  in  the  wing  spanwise  direction  has  minimal  effect  on  the  solution 
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(c)  d  =  0.700  (d)  d  =  0.850 

'    Figure  5.11  Cp  distribution  plot  (Grid  3) 
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accuracy.  It  is  also  noted  that  the  maximum  time  step  allowed  by  the  flow  solver 
is  about  At  =  0.20. 

The  final  grid  system,  Grid  set  5,  considered  in  this  study  is  regenerated 
from  Grid  set  3  by  halving  the  streamwise  wing-grid  spacings  and  the  fuselage 
cross-sectional  grid  spacings;  also,  five  additional  grid  lines  have  been  added  right 
next  to  the  wing  root.  The  refined  six  block  grids  obtained  again  have  no 
negative  Jacobian  when  it  is  evaluated  by  finite  difference  approximations.  In  the 
flow  simulation  process,  however,  the  maximum  time  step  allowed  by  the  flow 
solver  is  now  reduced  to  about  0.05,  which  makes  it  very  expensive  in  the 
simulation.  The  solution  process  has  been  carried  out  for  14,500  time  steps  which 
give  nondimensionalized  time  of  only  t^o,  =  672.  In  any  case,  the  predicted  Cp 
distribution  and  the  shock  wave  location  and  sharpness  on  the  lower  wing  surface 
are  in  better  agreement  with  measured  data,  as  shown  in  Figure  5.13,  than  the 
other  cases.  Again,  the  grid  refinement  has  not  improved  the  predicted  result  on 
the  upper  wing  surface. 

5.3  Complex  Flow  Physics 

In  this  section  the  complex  flow  phenomena  and  physics  of  the  flow 
problem  simulated  by  the  Reynolds-averaged  thin-layer  Navier-Stokes  equations 
with  the  use  of  Grid  set  5  are  presented  and  discussed  for  references.  As  shown 
in  Figure  5.13,  the  predicted  result  on  the  upper  wing  surface  is  not  quantitatively 
accurate;  however,  the  predicted  Cp  distribution  and  shock  wave  location  on  the 
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Figure  5.13  Cp  distribution  plot  (Grid  5) 
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lower  wing  surface  at  out  board  wing  cross-sections  are  in  good  agreement  with 
the  measured  data.  Further  insight  on  the  flow  characteristics  can  be  obtained 
from  Mach  contour  plots.  Figures  5.14  and  5.15  present  the  Mach  contour  plot 
with  AM  =  0.1  on  surfaces  a  short  averaging  distance  (Sg^)  from  the  wing 
surfaces.  The  simulated  results  show  that  the  flow  field  above  the  upper  wing 
surface  is  subsonic.  For  the  lower  wing  surface  flow  field,  the  location,  extent, 
and  strength  of  the  shock  wave  are  clearly  identified  by  the  Mach  contour  plots 
shown  in  Figure  5.15. 

The  shear  stress  on  the  wing  surfaces  also  has  been  computed  for  further 
understanding  on  the  flow  physics.  The  dashed  and  sohd  lines  in  Figure  5.16 
represent  the  streamwise  shear  stress  distribution  (t^j)  on  the  lower  and  upper 
wing  surfaces  at  four  wing  cross-sections.  On  the  lower  wing  surface,  a  sudden 

ft      .  T  :  ■ 

decrease  in  the  shear  stress  implies  a  sharp  increase  in  the  boundary  layer 
thickness;  also,  it  identifies  the  location  and  sharpness  of  the  shock  wave.  These 
shock  wave  related  results  are  in  good  agreement  with  those  shown  in  Figure  5.13 
and  Figure  5.15.  Moreover,  the  shear  stress  distribution  along  the  cross-section 
d  =  0.225  shows  that  there  is  a  very  small  shock  induced  separation  region. 


89 


(b)        =  0.68  (L  =  20) 
Figure  5.14  Mach  contour  plot  above  upper  wing  surface  (Grid  5) 
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(b)  S,^  =  0.68  (L  =  20) 
Figure  5.15  Mach  contour  plot  below  lower  wing  surface  (Grid  5) 
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(c)  d  =  0.700  (d)  d  =  0.850 

Figure  5.16  Streamwise  shear  stress  distribution  plot  (Grid  5) 


CHAPTER  VI 
SUMMARY 


A  multi-block  method  being  explored  is  further  developed  and  investigated 
for  the  simulation  of  Mach  0.8  transonic  turbulent  flow  past  a  wing-fuselage 
configuration  at  -3.0  degree  angle  of  attack.  In  this  method  the  flow  field  of 
interest  is  first  divided  into  six  contiguous  blocks  (one  nose  block,  two  side  blocks, 
two  wing  blocks,  and  one  tail  block)  in  such  a  way  that  each  block  is  partially 
bound  by  a  solid  surface.  Accordingly,  the  solid  surface  of  each  block  is  mapped 
onto  a  complete  rectangular  boundary  plane  in  the  computational  domain  for 
effective  calculation  of  eddy  viscosity  as  well  as  for  direct  application  of  a  thin- 
layer  Navier-Stokes  flow  solver.  Then,  in  the  solution  process,  each  block  is 
treated  as  an  independent  flow  problem,  modeled  by  the  unsteady  Reynolds- 
averaged  thin-layer  Navier-Stokes  equations,  with  interface  boundary  conditions 
updated  at  every  time  step.  The  turbulence  closure  model  used  is  the  Baldwin- 
Lomax  algebraic  eddy  viscosity  model. 

The  multi-block  method  being  investigated  has  a  distinct  advantage  in 
design  application  in  the  sense  that  a  local  change  to  the  configuration  requires 
only  the  related  block  grid  to  be  regenerated.  In  apphcation  the  method  is  indeed 
very  simple  and  straightforward;  however,  the  special  features  of  the  method 
have  resulted  in  excessive  distortion  on  block  domain  transformation  and 
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consequently  make  it  difficult  to  generate  good  block  grids.  Hence,  special 
measures  and  proper  techniques  for  block  grid  generation  are  investigated  and 
developed  in  this  study. 

A  one-dimensional  adaptive  grid  generation  technique  based  on  a 
variational  principle  has  been  employed  to  generate  wing  surfaces  adaptive  to 
measured  and  interpolated  pressure  gradients.  With  the  solid  surface  and  the 
corresponding  outer  surface  generated,  the  six  block  grids  are  generated  by  an 
algebraic  Two-Boundary  Method  based  on  transfinite  interpolation.  The  blending 
function  used  is  a  Hermite  cubic  polynomial  which  can  provide  good  grid 
characteristics  near  the  surface.  Also,  a  hyperbolic  tangent  clustering  function  has 
been  used  for  effective  grid  point  distribution  and  sufficient  grid  resolution  in  the 
viscous  boundary  layer. 

To  gain  further  insight  and  understanding  on  the  multi-block  method,  a 
base  set  of  adaptive  block  grids  having  the  dimension  of  55x25x60,  79x13x60, 
79x35x60,  79x35x60,  79x13x60,  and  51x25x60,  respectively,  for  Block  1  to  Block  6 
has  been  generated  for  use  in  the  flow  simulation.  These  block  grids  have  better 
grid  characteristics  in  the  sense  that  the  Jacobian  of  transformation  evaluated  by 
finite-difference  approximations  is  positive  at  every  grid  point. 

The  CPU-expensive  numerical  experiments  have  been  conducted  to 
investigate  the  effect  of  grid  resolution  to  the  solution  algorithm  and  accuracy. 
Important  findings  obtained  are: 
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(1)  Doubling  the  wing  spanwise  grid  resolution,  from  35  points  to  69  points 
has  little  improvement  on  the  predicted  wing  surface  pressure  distribution. 

(2)  Increasing  the  wing  chordwise  grid  points  from  79  to  157  does  improve  the 
solution  accuracy;  however,  the  maximum  allowable  time  step  is  drastically 
reduced  to  0.05. 

(3)  The  predicted  results  on  the  upper  wing  surface  agree  only  qualitatively 
with  the  measured  data,  which  could  be  due  to  poor  choice  of  the  outer 
boundary  location. 

It  appears  that  aided  by  the  present  multi-block  method,  if  appropriate 
number  of  grid  points  is  used,  the  predicted  pressure  distribution  and  the  shock 
wave  location  on  the  lower  wing  surface  can  be  in  good  agreement  with  the 
measured  data.  Based  on  this  observation,  it  is  concluded  that  the  multi-block 
method  is  a  very  promising  approach  which  can  be  further  developed  for  complex 
configuration  aerodynamics  simulation. 
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